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EXECUTIVE SUMMARY 

Ocean Data ^ ‘^tems, Inc. has performed advanced development of a set of 
atmospheric analysis and prediction models in support of the SEASAT Program of 
NASA under Contract NAS5- 24468 to Goddard Space Flight Center. 

Analysis Model Development 

The task objectives were to devise, implement and test modifications (to 
existing NASA-ODSl objective analysis models) to; (1) facilitate the use, 
distribution and coupling of ocean-area surface winds in the analysis models; and 
(2) examine satellite data dependencies in the analysis models. 

The analysis models, which utilize a 125x125 polar stereographic grid of 
the Northern Hemisphere, have been modified in order to incorporate and assess 
the impact of (real nr simulated) satellite data in the analysis of a two-day 
meteorological scenario in January 1979. Such program /procedural changes 
include: 

• a provision to utilize winds in the sea level pressure and multi-level 
height analyses (1000- IOC MBS); 

• the capability to perform a pre-analysis at two "control levels" (1000 
MBS and 250 MBS); 

• a greater degree of wind- and mass-field coupling , especially at these 
control levels; 

• an improved facility to bogus the analyses based on results of the pre- 
analysis; and 

• a provision to utilize (SIRS) satellite thickness values and cloud-motion 
vectors in the multi-!evel height analysis. 
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Once modified, the analysis models were tested (and adjusted further) using 
observations for January 1979 to establish workable procedures and parameter 
settings. Only the 19-21 January period is highlighted in this Final Report. 

Two types of experiments were conducted (and documented) to determine 
the potential impact of real (or simulated) satellite data on this set of models: 

• an ocean-area surface wind experiment to determine the extent to 
which SASS winds might be used to enhance the specification of sea- 
level pressure (and low-level flows, generally) during a two-day 
analysis sequence, and 

• an upper-level satellite information experiment to assess the impact 
oi cloud-motion vectors and SIRS soundings during a two-day analysis 
sequence. 

In the surface wind experiment, suitably-adjusted ship winds (used in lieu of 
unavailable SASS winds) were used, while associated ship pressures were with- 
held , in specified octants of the hemisphere. Location of octants (for with- 
holding pressure reports) depended on the analysis time. This analysis sequence 
ran from OOOOZ 19-21 January and led to Chart Set C. A second sequence for 
the same period in which both wind and pressure reports were withheld in the 
same time-dependent octants (Chart Set D) was executed to provide an upper 
bound on analysis error. The same upper-air data base was used in both 
sequences. Any differences between corresponding charts in Sets A, C and/or D 
are attributable to changes in the sea-level pressure analyses (and the set of 
twelve-hour forecasts used to bootstrap each analysis sequence). Chart Set A is 
the baseline analysis sequence which used ^ available observations. 

Results from such tests indicate that the surface winds had a positive 
impact on the SLP analyses, and reduced the potential analysis error by about 
25% in the two-day sequence. Such effects are propagated vertically into the 
entire troposphere through the distribution and coupling mechanisms in the 
models. The differences between Set C and Set A analyses tend to take on the 


spatial scales peculiar to the number and locations of withheld pressure reports 
(as would be expected). But, since SASS wind sets would be quite large 
(compared to the number and distribution of ship reports), we would expect that; 
(1) the differences would be of meteorological scale; and (2) the "potential" 
ar alysis error would be greatly reduced. 

Prediction Model Development 


The task objectives were; (1) to implement and test an iterative procedure 
(based on the Temperton scheme) in order to dynamically condition the initial 
state specification; (2) to modify that procedure and examine the effects on 
model behavior and program efficiency; (3) to modify the program to ingest new 
data at a second time origin and to dynamically condition new state (at the 
second time origin); and (4) to produce short-range forecasts to assess the impact 
of dynamic initialization. 

The (program) context for carrying out these task objectives was prepared 
and tested. A version of the five-layer, coarse-mesh (63x63) prediction model 
(PECHCV) was modified to permit the dynamic conditioning/adjustment [using an 
iterative procedure devised by Temperton (1973) and tested in a multi-layered 
model by Kesel and Wellck (1975)] of any state specification (i.e., the specifica- 
tion of model variables at any time origin). The new version, referred to as 
DYNMOD, has been coded to permit the selection (via data cards) of any number 
of iterations ("orbits") of any orbital length. For this study, the length consists 
of five (time) steps forward and five steps backward (about some time origin) per 
orbit. Even the "static initialization" case (zero orbits) can be handled by 
DYNMOD, since aU of the dynamic initialization runs begin with static initiali- 
zation. 

Analysis model outputs (for 19/OOOOZ, 19/1200Z, and 20/0000Z January 
1979) were used for the static initialization phase of the forecasts. In some 
forecasts this static initialization was augmented by a specified number of 
Temperton-style orbits before commencing the two-day forecasts (dynamic, 
initialization at a single time origin). In still other forecasts, grldded values (for 


a pie-shaped region in the Pacific) from an analysis of conditions twelve hours 
into the forecast period were ingested/assimilated and dynamically adjusted 
using additional orbits before completing the remainder of the two-day forecast 
(i.e., dynamic conditioning at two time origins with data assimilation at the 
second). Identifiers for these forecasts are as follows: 


Run 

Number of 

Initial Conditions for 

Type 

Orbits 

19/OOOOZ 

19/1200Z 

20/ OOOZ 

static 

initialization 

only 

0 

E 

- 

- 

dynamic 
conditioning 
at a single 

6 

1 

F 

1 

- 

time origin 

12 

G1 

G2 

G3 

dynamic 
conditioning 
at two time 

12+12=24 

HI 

H2 

! 

H3 

origins plus data 
assimilation 



i 



All forecasts were compared to the appropriate analyses within the "baseline" 
analysis sequence (which used all available observations). Additionally, model 
behavior was monitored via a set of diagnostics computed each forecast hour. 
The most useful measure of dynamic conditioning of model variables is the 
(squared) mass divergence (at each model level). 


The results from these forecasts are described in terms of program 
efficiency, model behavior and forecast quality. 


• Program Efficiency 

The FTN (Fortran) program DYrN'MOD is modular, quite general, and 
substantially transportable. It was executed on both a CDC CYBER 175 (using 
NOS/BE) and a CYBE.v m (using SCOPE). Central memory (user) requirement is 




sac/'.^aj! 


65K (decimal) words. No Extended Core Storage (ECS) was used. About three 
million words of mass storage are needed for each two-day forecast. 

Dynamic initialization may be performed using any number of orbits 
(including zero orbits) of any orbit length. Data assimilation in a specified 
geographical region and dynamic conditioning at a second time origin are 
"hardwired" in the code at present. With only minor changes thereto, other 
regions and/or other time origins could be handled. 

On a CYBER 74, static initialization takes 3,770 CP-seconds and 7,460 
seconds of I/O. There is a substantial overlap of CP and I/O during model 
integrations, but not during initialization and output preparation. Each orbit 
adds about 190 CP-seconds and 400 I/O-seconds to these totals. 

• Model Behavior 


Dynamic adjustment/conditioning of any state specification does not 
enhance the validity of that specification. Rather, it enhances the quality of the 
relationships between the mass fields and winds, especially in the numerical 
sense. Through numerical/dynamic consistency, initialization shock is reduced 
and modeling of secondary physical effects is improved (in the first several 
forecast hours). 

Examination of the (squared) mass divergence statistic in Forecasts E, 
F, and G1 shows that: (1) the reduction of initialization shock is proportional to 
the number of orbits made; (2) dynamic initialization greatly accelerates the 
geostrophic adjustment process in the model — a task which would take the 
model itself from one to two forecast days to accomplish with its own 
"adjustive" powers; (3) the use of twelve orbits was satisfactory for this study. 

Dynamic conditioning appears to have less impact on the layer-mean 
kinetic energy (in these two-day forecasts). The differences noted between the 
static and dynamic (initialization) cases tend to be small (at the outset) and 
remain small throughout the forecasts. The adjustments noted tended to be in 
the right direction (i.e., closer to values the forecast ultimately produces). 


There is an impact on the model's precipitation mechanism. Convec- 
tive precipitation tends to be reduced (i.e., the size of the affected area and the 
average amounts). Large-scale precipitation appears to be increased by about 
10-20 percent. (Vertical velocity patterns become more mature during the 
adjustment process — giving rise to more precipitation in the first forecast day.) 

• Forecast Quality 

Differences in the type/complexity of the initialization did not lead to 
large differences in the two-day forecasts which were produced. This appears to 
be true in spite of the nurnerical/physical benefits juct discussed. 

Forecast Runs E, F and G1 (which all used initial conditions for 00002 
19 January) produced RMSE scores at 500 MBS of 46.2 - 46.4 meters, and 4.97 - 
4.99 MBS at sea level. Forecast HI, which ingested new (analysis) data for 
19/12002, had 46.4 meters and 5.01 MBS at 500 MBS and sea level, respectively. 
Data ingestion did not improve the forecast error. 

Forecasts G2 and H2, initialized at 19/12002, differed in RMSE scores 
by only 0.05 MBS and 0.9 meters at sea level and 500 MBS, respectively. Thus, 
data assimilation produced a positive, but small impact on the forecast. 

Forecasts G3 and H3, initialized at 20/00002, differed in RMSE scores 
by only 0.06 MBS and 1.4 meters — with data assimilation producing a positive 
but small impact. 

An examination of the central pressures and pressure profiles of major 
low-pressure systems (handled by these forecasts) indicates a minor (and mixed) 
impact when using dynamic conditioning and/or data assimilation. Some features 
benefitted and others did not benefit. 

We have tentatively concluded that the computational expense is 
difficult to justify under the conditions of this study. The lack of precise 
knowledge of initial conditions over the world's oceans is a much larger problem 
than problems related to numerical consistency. 
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I. 


INTRODUCTION 


Ocean Data Systems, Inc. (ODSl) has designed, developed and tested 
atmospheric analysis and prediction models of varying (grid) resolution in support 
of the SEASAT Program of the National Aeronautics and Space Administration 
(NASA); (I) under Contract NASW-2558 to Econ, Inc.; and (2) ur.Jci Contract 
No. 954668 to the 3et Propulsion Laboratory OPL). Some advanced development 
of these models was performed under Contract No. NAS5-24469 to Goddard 
Space Flight Center (GSFC). In the latter effort, procedures for enhancing the 
computational viability of SEASAT-type data were identified and tested. j|.>®ci- 
fically, the effects of three discretionary procedures in objective analysis were 
studied: (1) varying the weights in the Patvern Conserving Technique; (2) varying 
the influence functions for observations: and (3) examining the effects of using 
wind information in analyses of mass-structure variables. 

Finally, this Technical Report covers additional advanced development of 
the models under Contract No. NAS5-24468 to GSFC. The effort consisted of: 
(1) modifying the analysis programs and determining the impact of (real or 
simulated) satellite observations on the analysis of a meteorological scenario in 
January 1979; and (2) developing and testing a primitive-equation prediction 
model which is initialized using dynamic methods, and which is capable of both 
data assimilation and dynamic re-initialization at a second time origin. These 
advanced development efforts have focused on ways to prepare a better "model 
context" for assessing the utility of SEASAT-type observations in analysis and 
prediction. Ways were sought to ingest and to distribute such information into 
(spatial) scales which are less likely to get "washed out" in the adjustment period 
of a short-range prediction using a primitive-equation model. The objective 
analysis program sequence was modified rather extensively in order to facilitate 
the use of both satellite-sourced winds and temperature profiles, as well as the 
coupling of wind and mass-structure information. 

The testing of the analysis and prediction models was carried out with a 
real data set for January 1979 made available to ODSI by Fleet Numerical 
Weather Central, Monterey. For purposes of this Final Report, however, onlv 
the observations and computer runs for the period 19-22 January were used. 
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FNWC also provided some computer support in the early phases of this 
contractual effort under an agreement between FNWC and NASA. As the work 
progressed to those phases requiring more frequent turnaround and longer 
production runs, it became necessary to shift the data base and computer 
programs to other computational facilities. This led to delays and to additional 
work. Nevertheless, the work has been successfully completed. This transition 
to other computer systems was facilitated, in part, by FNWC which provided to 
Onsi (through NASA) some of the system software peculiar to the FNWC system 
environment. 


Section II pertains to analysis model development. There is a description 
of the models, and the descriptions and test results of a marine-wind experiment 
and an upper-level satellite-information experiment. Section III provides the 
descriptions of and test results from use of prediction models initialized using 
conventional (static) and dynamic methods. Comparative results for one- and 
two-time-origin dynamic initialization model executions are provided. Sections 
IV, V and VI contain the charts, tables and figures, respectively. Five appendices 
are also provided. Appendices A and B provide a ready reference for readers 
interested in the use of PCT in scalar and vector analyses, -espectively. The 
derivation of the Mass- Structure Linear Transformations is provided in Appendix 
C. Appendix D describes the procedure used to merge/check mandatory and 
significant level radiosonde data. A brief description of the procedure used to 
obtain temperatures and heights above the 100 MB level is contained in Appendix 




II. ANALYSIS MODEL DEVELOPMENT 


The objective analysis models developed by ODSI for the SEASAT studies 
employ the Pattern Conserving Technique (PCT). The models and procedures 
were modified rather extensively in order to accomplish the proposed objectives 
in this effort. The description of PCT is provided in Part A, as well as in 
Appendices A and B to this Report. Part A also: (1) outlines the program 
sequence used to complete the several analysis production tests; and (2), 
describes the modifications to models and procedures necessary to accomplish 
study objectives. Part B describes the impact of real (and/or simulated) satellite 
observations on the analyses from a two-day scenario in January 1979. Part B is 
presented in three parts. Part 1 introduces the "baseline" analysis sequence 
produced using all available data. Part 2 describes a marine-wind experiment in 
which shipboard surface winds (used in lieu of SASS winds) are employed to 
determine the extent to which ocean-area surface winds may be used to help 
maintain the sea-level pressure specification in a two-day sequence. Part 3 
describes an upper-level satellite information experiment in which the impact of 
cloud-motion vectors and SIRS reports is determined (in the same scenario). 

A. Description of Models 

1. The Pattern Conserving Technique 

The procedure/method used by ODSI for all analysis modeling in 
support of SEASAT is called the Pattern Conserving Technique (PCT). It has 
been used by ODSI to analyze sea-surface temperature; sea level pressure; and 
the multi-level temperatures, heights and winds (from 1000-100 MBS). It is a 
multi-cycle procedure in which each cycle consists of: 

• Assembling the observations to grid points. 

• Solving the PCT minimization equation. 

• Re-evaluating the weight of each observation. 

• Adjusting the influence function for each observation. 


The method lends itself very well to the "engineering adjust- 
ments" which have to be accomplished in operational analysis. With an 
appropriate selection of weights, control can be exercised over which of the 
characteristics will be emphasized in the final analysis. In the assembly stage, 
one exercises control over the spatial scales being affected by each datum. The 
relative importance of each observation — as it changes in each cycle — may 
also be controlled. The complexity of the minimization equation may be tailored 
to the purpose of the analysis. (In this contract effort, for example, the pressure 
and height analyses were modified to use wind information.) Filtering may be 
adjusted as desired, .,.'en though the "perfect filter" does not exist. 

A rather complete description of the models (the PCT method; 
relevant equations; computer program descriptions) is available in Volume II of 
the ODSI Final Report to the 3et Propulsion Laboratory (dated 30 September 
1977) under Contract 954668. Additional model discussion/description is con- 
tained in the ODSI Final Report to Goddard Space Flight Center (dated 
November 1978) under Contract NAS5-24469. Appendices A and B to this report 
explain the use of PCT in the analysis of scalars and vectors, respectively, for 
readers without access to the aforementioned references. 

2. Analysis ProRram Sequence 

All of the analysis results presented in this report were generated 
using a 125x125 grid-point polar stereographic projection. This grid has a 
resolution of 190.5 kilometers at 60*^ latitude. This resolution should be better 
able to accommodate the spatial scales observed by SEASAT, yet is not as 
computationally burdensome as the 187x187 grid used by ODSI previously. 

Figure II- 1 depicts the order of execution of the jobs required to 
complete one analysis sequence. The jobs are executed serially since the output 
from one program is often required input for the next program in the sequence. 
The sequence is initiated every 12 hours with the analysis date and hour specified 
as the current "date-time group" (DTG). Once the DTG is set and the required 
data records are made available, a first guess field for a -6 hour (6 hours before 
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the current DTG) sea level pressure (SLP) analysis is generated. To prevent time 
lag, this guess field is time-centered, combining the 12 hour old SLP analysis and 
the 12 hour SLP forecast from the last sequence run. With this guess field, the - 
6 hour SLP analysis is performed using the PCT as described in Appendix A. The 
resulting field serves as the major input for the ontime SLP guess field which is 
produced in the next program. At the same time, the guess fields are generated 
for the preliminary upper air temperature ana'ysis. The time-centered tempera- 
ture guess fields represent a blend of the 12 hour old temperature analyses and 
the previous run's 24 hour forecast temperature fields. While the guess program 
is executing, the sea surface temperature (SST) analysis can also be run since its 
execution is independent of other analyses. With the completion of the guess 
field program, the preliminary ontime SLP analysis is initiated. This is followed 
by a program which checks the validity of the radiosonde data. (This radiosonde 
checking program is described in Appendix D). Next, the preliminary upper air 
analysis is performed. [The preliminary upper air analysis produces upper air 
temperature and height guess fields for the final analysis which are consistent 
with two "control" levels. The 1000 MB and 250 MB levels were selected because 
they generally have the best data coverage. The 1000 MB level is derived from 
the data-rich sea level pressure analysis, and 250 MB is the level where both 
aircraft reports and satellite cloud-motion vectors are conce.itrated. Deriving 
the other guess fields from the analyses at these two levels provides a consistent 
and vertically coupled guess field structure.] There is now a complete set of 
guess f'elds and the final combined analysis (final SLP and upper air tempera- 
ture, height, and wind analyses) can be performed. To review, the key steps in 
the upper air analyses, preliminary and final, are: 

• The 1000 MB height guess is generated using the "prelimi- 
nary" SLP analysis and the preliminary 1000 MB temperature 
guess. 

• The 1000 MB height analysis is performed. 

• The preliminary first-guess height fields (8) for the 950-250 
MB levels are generated using the 1000 MB height and a 
hydrostatic integration of the preliminary first-guess temper- 
atures. 
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• The preliminary 250 MB height analysis is performed, using 
the bogus capability as I'ecessary. 

• The preliminary D-values for the seven (7) pressure levels 
between 1000 and 250 MBS are computed by jsng a In p 
interpolation of the preliminary D-values at 1000 and 250 
MBS. 

• The twelve (12) preliminary heights (1000-100 MBS) are then 
"transformed” into ten (10) static stability values — which 
are then adjusted as necessary. 

I^Note: The reader should see Appendix C for details of the Mass Structure 

Linear Transformations being used. In order to retrieve twelve heights (or 
temperatures) from these static stabilities, two additional conditions must be 
specified: we have chosen the 1000 MB height and the 1000-250 MB thickness.^ 

• The twelve (12) first-guess temperature fields are retrieved 
from the static stabilities, 1000 MB height and 1000-250 MB 
thicknesses. 

• The final SLP analysis is performed, bogusing as necessary. 

• The upper-air temperature analysis is completed using guess 
fields just retrieved. 

• A 1000 MB first-guess is generated from the final SLP and 
the 1000 MB temperature fields. 

• The height analysis is performed one level/layer at a time, 
using the heights at the level just analyzed and the final 
analyzed temperatures (for the given layer) in a hydrostatic 
integration to produce the guess heights for the next level to 
be analyzed. 
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• Usinj* regression equations described in Appendix E, the 
temperature/height values are extended above 100 MBS (to 
10 MBS). 

• The winds at twelve (12) levels are analyzed using the PCT 
method described In Appendix B. 

A five-level 63x63 grid-point primitive equation forecast model is 
used to bootstrap the analysis sequences. To initialize the forecast model, the 
required 125x125 analysis fields must be "dezoomed” to 63x63 fields. The model 
(without dynamic initialization) is integrated to 24 hours to produce the required 
forecast input fields for the next analysis. Since forecast velocity divergence is 
required as a constraint in the wind analysis, a separate program is -un to create 
these fields. 


3. Program/Procedural Modifications 

A considerable amount of u;i>e and thought was given to possible 
approaches to the proposed analysis model modifications. The design which was 
selected provided the following essential features; (1) it was capable of utilizing 
all available satellite (and conventional) observations; (2) it emphasized the use 
of ocean-area surface winds, of cloud-motion vectors in the 300-200 MB layer, 
and of satellite temperature profiles; (3) it was capable of using wind informa- 
tion in the sea-level pressure and multi-level height analyses; (4) its complexity 
was numerically and dynamically consistent with the requirements and charac- 
teristics of the prediction model; (5) it was computationally economical; (6) the 
work could be successfully accomplished within project resource limitations; (7) 
it had the potential for fulfilling the task objectives; i.e., to distribute the 
reported information and to couple the estimates of the winds and mass- 
structure parameters. 

Several types of analysis program/procedural modifications were 
incorporated for this contractual effort. Brief descriptions are listed below. In 
subsequent paragraphs, each modification will be described more fully. 
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• analysis model re-design allowing for a pre-analysis at the 
1000 MB and 250 MB control levels — with a facility for 
bogusing as necessary. 

• provision to utilize wind informatiw (gradients) in the sea- 
level pressure and multi-level height analyses. 

• changes in the (data) assembly procedure and in the number 
of analysis cycles (in the sea-level pressure and multi-level 
temperature ^alyses) to better accommodate and distribute 
observed information. 

• provision to return (functional values) to long-term clima- 
tology in regions of few or no recent observations — in order 
to control the analyses. 

• time-centered blending of first-guess fields, using a recent- 
past analysis and a short-range (PE) forecast. 

• changes in the type and degree of filtering (of model out- 
puts). 

The most significant modification of the analysis program set 
involved a major re-design which permits a pre-analysis at the 1000 MB and 250 
MB "control levels" (new feature). At the control levels, which tend to be data 
rich anyway, use of the bogusing capability r» encouraged and facilitated in order 
to arrive at the best possible specification of analysis variables. The thickness 
values from satellite (height) observations are relied upon to produce the 
(important) 250 MB height pre-analysis. The 1000 MB height first-guess field is 
generated from a (very) data-rich sea-level pressure analysis. Very often, the 
final 1000 MB heights are subject to change with the addition of 500-600 
radiosondes and with the (different) procedures that are embedded in the upper- 
air analysis programs. Winds are allowed to strongly influence the height 
patterns at both control levels. The D-values produced at the control levels in 
the pre-analysis are then distributed to remaining (intermediate) pressure levels 
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using In p interpolation. C<y>sistent temperatures are generated from this height 
structure using the linear transformatiwis described In Appendix C. 

Second, the PCT approach proved to be very useful when it was 
decided that wind information had to be incorporated in pressure and (multi- 
level) height analyses. Each wind report Is checked against a gross tolerance 
which is similar to that being used in the wind analysis. If accepted, the 
velocities are then converted to height/pressure (gradients. These gradients are 
substituted for the (height/pressure) guess-field gradients at the nine grid points 
(3x3) closest to the wind observation. In this manner, the gradients act as a 
constraint on the solution of the PCT minimization equation(s). ( ine calcula- 
tions are made in the subroutine BKGRND.) The (3x3) stencil of grid points 
mentioned above could not be enlarged without a major recoding effort because 
of the central memory limitation. Otherwise, the spatial scales being influenced 
in this way might have been quite different. 

The assembly procedure (used to get the observed values to the 
grid points) was also modified: (1) to better distribute observed information; and 
(2/ to make the procedure more efficient and still meet the computer resource 
limitations. For the most part, the procedure which was used (previously) with 
the 63x63 models still applies. The exceptions are noted below: 

• Observations south of 5°N are not used (they cause more 
problems than they solve). 

• For the 125x125 grid, observations are sorted according to 

the "3-row" (3 = 1,1 ), and then grouped into five-row sets for 

inclusion/assembly. The central memory restriction dictates 
that only one such five-row set can be assembled at a time. 

• The region influenced by each observation (now) depends 
upon; data density; cycle number; analysis type; and (in the 
SLP and S5T) the differential properties of the guess field. 
Generally, a datum doesn't affect grid points very much if 
such points are already affected by closer and/or more 
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oensely-spaced observations. The influence region gets 
smaller with each cycle through the analysis. The region 
being affected by a report is larger in the upper-air analyses 
than at sea -level, other factors being the same. Regions with 
strong gradients and/or curvature cause the influence radius 
to be reduced according to the magnitudes of such gradients 
and/or Laplacians, but the provision applies only to the 
analysis of sea-level pressure and sea-surface temperature. 
Finally, the toss-out criterion also depends (now) on latitude - 
with the tolerance being reduced as the latitude decreases. 

• Each analysis (except for the wind analysis) now returns to 
long-term climatology in regions for which there are few or 
no recent observations. North of 20°N, the specification 
returns at a 5% rate each analysis time. Below ?0°N the rate 
of return is invei-sely proportional to latitude. The informa- 
tion density, of course, is the primary consideration in the 
calculation. 

The analysis programs have been modified so that they may be re- 
starte d from a given DTG ~ but with slightly altered procedures and constants 
for the first such DTG in the new analysis sequence. Such alterations consist of: 

• an increase in the data reject limits. 

• extra smoothing and filtering of outputs (since the first guess 
for the re-start is climatology). 

• altered PCT weights in the sea-level pressure and height 
analyses (these analyses use wind gradients which enter via 
the PCT equations). 

• data weights are not re-evaluated after each analysis cycle. 

• the PCT stage is by-passed in the U/A temperature analyses. 
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G. impact of Satellite Observations on Analysis Structures 


After describing the ''baseline" analysis sequence, the remaining 
materia! is presented in two major parts: (1) an ocean-isiea surface wind 

experiment; and (2) an upper- level satellite information experiment. Each type 
of result may be compared to that obtained in the baseline analysis sequence. 

In the part of the study involving objective analysis, there were two 
proposed objectives: (1) to implement and test an analysis procedure to 

distribute surface information vertically into the lower troposphere, and to 
couple the winds and mass-structure variables; and (2) to examine satellite data 
dependencies in the analysis models. These experiments were conducted in order 
to demonstrate the proposed objectives. 

Four separate two-day jurialysis sequences were produced. Each 
sequence is presented with its own chart set. The key features of each analysis 
chart set are defined below; 

Set A: Baseline analyses made using all available data. 

Set B: Analyses made without cloud-motion vectors and SIRS 

reports. 

Set C: Analyses in which shipboard pressure reports are withheld (in 

specified octants) but associated ship winds are used. 

Set D: Analyses in which entire ship reports are withheld (in speci- 

fied octants), 

I. Baseline Analyses 

The set of ODSl-NASA objective analysis models (programs) was 
used to produce analyses for a scenario extending from 00002, 19 January, 
through OOOOZ, 22 January 1979. The "baseline" set of analyses utilized all 
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available observations. This included; sea surface temperature; sea level 
pressure; and the temperatures, heights and winds at twelve constant-pressure 
levels from 1000 to 100 i\tBS. The heights and temperatures at 50 MBS, used to 
initialize the prediction model, are produced using regression equations described 
in Appendix E. Winds are derived from the heights at the 50 MB level. 

Although Versatec charts were produced for most analysis param- 
eters and levels, only a portion of them is included in this report. Charts IV- 1 
through IV-26, referred to as Chart Set A, contain the relevant information in 
the baseline analysis set. These should be examined by the reader in order to 
understand the meteorological context in which the experiments w^?re conducted. 

Charts lV-20 through IV-22 contain the 48-hour analysis change 
patterns at sea level and 500 MBS. In this scenario one can see some very large 
(and interesting) changes occurring everywhere except in southeast Asia and 
central Europe. In IV-20, for example, there are many change centers with 
magnitudes greater than thirty (30) millibars in this two-day period. The chart 
discussion will tend to focus on the Pacific area. Thus, the reader should 
examine the evolution of pressure systems shown in Charts IV-1, IV-8, IV- 10, and 
IV- 12, in that region. 

In Section II, the emphasis is placed on sea level pressure, and the 
temperature and wind at both 900 MBS and 250 MBS, levels at which cloud- 
motion vectors seem to be most plentiful. In Section III, the emphasis is placed 
on sea level pressure and 500 MB height, levels at which prediction models are 
often evaluated. Thus, Chart Set A includes analyses which will be used to 
discuss the results of both vhe analysis and prediction tasks. 

2. Marine-Wind Experiment 

The discussion is presented in two parts. First, there is a 
description of the experiment and the (additional) special modifications made in 
order to accomplish the proposed objectives. Then, the results of the analysis 
test demonstration are presented. 
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a. Description 


In this analysis task, the objective was to devise, implement 
and test appropriate modifications to the objective analysis models (under 
development for NASA by ODSI) so that SASS winds (and other satellite data) 
could be utilized to enhance the specification of the initial conditions (in ocean 
areas) for short-range atmospheric prediction. Specifically, ODSI proposed to 
devise and test a method for distributing surface-level information vertically, 
into the lower troposphere and for coupling the available (estimates of) wind and 
mass-structure parameters. In addition, ODSI proposed to examine what may be 
termed "data dependencies" in this particular set of analysis programs/models. 

Having referred (earlier) to limitations on project resources, 
it is important to note that the initially-envisioned comprehensive data set 
(containing SASS winds) did not materialize. Thus, ODSI decided to simulate 
SASS winds with ship (surface) winds — and to design experiment to determine 
the extent to which ocean-area surface winds may be used to maintain/enhance 
the sea-level pressure specification (as well as the flows in the lower tropo- 
sphere). Ship winds and pressures had to be withheld (say) from certain octants 
to define an "upper bound" on analysis error in one two-day analysis sequence; 
but only ship pressures would be withheld from the same octants in a second 
analysis sequence for the same two-day period. The octant(s) in which the ship 
reports were handled in this manner represented a simulation of the time-varying 
areal coverage of satellite data. Finally, the ship winds had to be "adjusted" in 
an appropriate way to make them more acceptable to the experiment. 


In the sea-level pressure analysis, for example, the geo- 
strophic wind law was altered in order to simulate the reduction of speed due to 
friction and the turning of the wind across the isobars. The coefficients, R and 
S, control the speed and direction, respectively, in the final version of the wind 
law shown below: 


u = 




V = 


1 

M 



S 
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where M = Rp (f + S ). Wind observations must first pass a gross-error check in 
which they are compared to winds derived from the first-guess pressure 
distribution using [ll.l] and [11.2]. If the speed difference exceeds a specified 
fraction of the guess-field speed, the observation is not used. If it is accepted, 
the pressure gradients at the nearest grid point(s) are set to the values implied 
by that observation. If more than one observation influences a grid point, the 
implied values from the (several) reports are averaged. With appropriate tuning 
of the PCT equation constants, an optimum mix of assembled (functional) values 
and differential properties may be obtained. 

b. Results of Tests 

Chart Sets C and D will be referred to in this discussion. 
Chart Set C contains selected charts from the sequence of analyses in which ship 
winds were used — but associated ship pressures were withheld in specified 
octants. The objective was to simulate the extent to which SASS winds might be 
used to assist in the specification of the sea-level pressure. Through coupling 
mechanisms in the three-dimensional analysis, effects should also be discernible 
at and above 900 MBS. In Chart Set D, both the pressures and winds were 
withheld in the same specified octants. This analysis set may provide an 
upper bound to the analysis error in this scenario which runs from OOOOZ 19 
January through 00002 21 January 1979. Thus, there were nine (9) sea-level 
pressure analyses (spaced every six hours) and five (5) upper-level analyses 
(spaced every twelve hours) in each of the two sequences (Set C and Set D). 
Finally, it is also possible to assess the interactive/cumulative effects within 
each analysis sequence since each sequence is connected by its own set of (12- 
hour) PE forecasts which produce the necessary first-guess fields. 

Chart IV-47 (from Set C) contains the sea-level pressure 
analysis for 00002 19 January. About 280 (ship) pressure reports (6% of total) 
were withheld from this analysis. From these 280, there were 232 wind observa- 
tions accepted by the analysis. This analysis should be compared to the 
corresponding (Set A) baseline analysis provided as Chart IV- 1. The difference 








between these Set C and Set A analyses is provided as Chart lV-53. At OOOOZ, 
the octants being affected are in the eastern Atlantic and western Pacific. In 
the Atlantic the differences are generally smaller than three (3) millibars — but 
in the Pacific they tend to be much larger (8 millibars in one place). 

Charts lV-48 and lV-49 contain the (Set C) 900 MB Tempera- 
ture and Wind Analyses, respectively, for OOOOZ 19 January. Compare these to 
Charts IV-2 and IV-3 in the baseline sequence and examine the differences. 
Charts lV-54 and IV-55 actually contain the corresponding differences for the 
900 MB temperatures and winds, respectively. 

Next, we will discuss the charts two days into the (inter- 
active) analysis sequence. Chart IV-50 (Set C) contains the sea-level pressure 
analysis for OOOOZ 21 January. About 314 ship pressure reports (6.7% of total) 
were withheld from this analysis. Of these 314 reports, there were 236 wind 
observations accepted into the analysis. Compare this analysis to the corres- 
ponding analysis from the baseline analysis set shown in Chart IV- 12. The 
difference between the Set C and Set A analyses is shown in Chart IV- 56. 
Although it is difficult to determine visually, the (RMS) differences on 21 
January ure slightly smaller than on 19 January (0.55 vs. 0.62). Indeed, whatever 
cumulative differences there might be are not discernible in either the differ- 
ence charts or statistics. The differences at each analysis time are primarily 
related to the number/location of pressure values being withheld, as well as to 
the viability of spatial scales portrayed in each analysis. 

Charts 1V-5I and lV-52 (from Set C) contain the 900 MB 
Temperature and Wind Analyses, respectively, for OOOOZ 21 January. Compare 
these to the corresponding baseline analyses shown in Charts IV- 1 4 and IV- 15. 
The differences between corresponding 900 MB analyses are provided as Charts 
IV-57 and lV-58. Table V-I contains the RMS differences between Set C and Set 
A analyses in the sequence. The RMS temperature difference at 900 MBS varied 
from 0.02-0.04 degrees Celsius during the two-day sequence; and the correspond- 
ing RMS height difference at 900 MBS varied from 2. 1-3.4 meters. Also, the 
RMS wind component difference varied from 0.46-0.71 meters per second. Note 
that the differences extend to the 250 MB level, where they are about half as 
large as at 900 MBS, generally. 
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The charts in Set D may be used as a control — an upper 
bound on the analysis error in that analysis sequence from which entire ship 
reports (pressures and winds) were withheld in the same octants used in Set C. 

Chart IV- 59 (Set D) shows the sea-level pressure analysis for 
OOOOZ 19 January. Compare this to Chart lV-47 (from Set C) and to Chart IV-1 
(from Set A). Chart IV-65 contains the difference between the Set D and Set A 
sea level pressure analyses at this time. This may be compared to Chart IV-53 
containing the difference between Set C and Set A. Or, finally, examine Charts 
IV-71 and lV-72 which show the differences between Set D and Set C analyses on 
19 and 21 January, respectively. The spatial scale of these differences reflects 
the number and location of withheld reports, primarily. 

Charts IV-60 and lV-61 (Set D) show the 900 MB Temperature 
and Wind Analyses, respectively, for OOOOZ 19 January. Compare these to the 
corresponding charts in the baseline analyses (Charts IV-2 and IV-3), as well as to 
the corresponding charts in Set C (Charts IV-48 and lV-49). The differences 
between Set D and Set A analyses are contained in Charts IV-66 and IV-67. 

Next, we will introduce and discuss the (Set D) analyses two 
days into the (interactive) sequence. Chart IV-62 contains the sea-level pressure 
analysis for OOOOZ 21 January. The difference between this analysis and the 
baseline analysis is shown in Chart lV-68. Table V-2 contains the RMS 
differences between Set D and Set A analyses for several parameters/levels. 
[Generally, one would expect that the (Set D-Set A) differences would be greater 
than the (Set C-Set A) differences. A comparison of Table V-2 with Table V-1 
shows that such an expectation is realized in two-thirds of the parameters/leveis 
tabulated.] The sea-level pressure (RMS) values in the two tables indicate the 
use of the winds reduced the analysis error by about 25 per cent. The greater 
number of SASS winds, however, would lead us to expect a much greater impact 
on (operational) sea-level pressure analyses. Secondly, there is little indication 
of cumulative effects in these sequences. Had the (initial) differences been in 
larger spatial scales (as they would be with SASS winds), we would expect to see 
more cumulative (interactive) effects. 
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Charts lV-63 and IV-64 ccmtain the 900 MB Temperature and 
Wind Analyses, respectively, for OOOOZ 21 January. Compare these analyses to 
the corresponding charts in Set A (Charts IV- 14 and IV- 15) and in Set C (Charts 
lV-51 and IV- 52). The 900 MB temperature and wind differences (Set D minus 
Set A) are provided as Charts IV-69 and IV-70, respectively. 

3. Upper-Level Satellite Information Experiment 

As discussed previously, model modifications were made which set 
the stage for this sensitivity study employing the three-dimensional temperature, 
height and wind analysis codes. In this experiment, the main objective was to 
assess the impact of satellite temperature /height profiles cUid cloud-motion 
vectors on a two-day analysis sequence. 

a. Description 

The January 19-21 period was not only meteorologically 
interesting, but one in which the numbers of observations were greater than 
usual. Table V-3 shows that the numbers of radiosonde/rawinsonde soundings 
available for the upper-air analyses varied from 550 to 590, with an additional 50 
to 85 SIRS profiles. Table V-4 contains the numbers of cloud-motion vectors, 
which varied from about 110 to 260 with the typical concentrations both at 900 
MBS and in the 300-200 MB /one. (The level of each wind value is specified to 
the nearest ten millibars.) Table V-5 indicates the even larger numbers of 
aircraft wind reports. These were useful for proper specification of high-level 
flows, generally, but extremely so at the 250 MB control level. Although the 
impact of cloud-motion vectors was both measurable and positive (in this 
experiment), it would probably have been greater were it not for the fact that 
aircraft winds had already played such a large role in specifying flows at high 
levels. 

SIRS soundings presently contain heights from 850 MBS 
through 100 MBS. Since these heights are referenced to a NOAA control field, 
the absolute heights are of less value to these (ODSI) analyses than the implied 
thickness values. The ODSI 850 MB height, therefore, is employed to compute 
heights which are consistent. 
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b. Results of Tests 


Chart Sets A and B will be referred to in this discussion. 
Chart Set A contains selected charts from the baseline analysis sequence -- 
produced using ail available observations. Chart Set B contains corresponding 
analyses from the separate analysis sequence from which both SIRS soundings 
and satellite cloud-motion vectors were withheld. The two-day sequence began 
with observations for OOOOZ 19 January. During this sequence, the corresponding 
sea-level pressure analyses for Sets A and B will be more or less identical. But, 
the upper-air temperature, height and wind analyses for the two sequences will 
diverge. The following discussion will concentrate on the differences at 900 MBS 
and 250 MBS. 


Charts lV-27 and lV-28 contain the (Set B) 900 MB Tempera- 
ture and Wind Analyses, respectively, for OOOOZ 19 January. Compare these to 
Charts lV-2 and lV-3, respectively, in the baseline analysis sequence. The 
differences between Sets B and A are provided as Charts lV-38 and lV-39. On 
Chart lV-38, note three areas where the 900 MB temperatures differed by more 
than one degree Celsius. These patterns are of meteorological scale. There are 
two additional places where differences exceed 0.5 degrees Celsius. On Chart 
lV-39, the wind differences tend to be less than ten knots. 

Charts lV-32 and lV-33 contain the (Set B) 900 MB Tempera- 
ture and Wind Analyses, respectively, for OOOOZ 21 January. Compare these to 
Charts IV- 14 and IV- 15, respectively, in the baseline sequence. The Set B-Set A 
corresponding differences are provided as Charts IV-43 and IV-44. The magni- 
tudes of temperature difference patterns on lV-43 are in excess of two degrees 
Celsius in the central and eastern Pacific. The wind differences shown on IV-44, 
however, do not appear to have increased beyond those on 19 January discussed 
on Chart IV-39. 


Table V-6 contains the RMS differences between Set B and 
Set A for several analysis times and parameters. The RMS temperature 
difference at 900 MBS is 0.17 degrees on 19 January — and increases to 0.32 
degrees on 20 January. No further increase occurred on 21 January. (The same 
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result was observed at 250 MBS wherein the RMS temperature difference 
doubled (to 0.22 degrees) after one day, but did not increase much beyond that 
amount in the second day.) Table V-6 also shows that the RMS wind differences 
at 900 MBS did not increase beyond 0.6 m/sec and 0.3 m/sec for the u-and v- 
components, respectively, as shown on the first day. 

Chart lV-29 contains the 250 MB Height Analysis for OOOOZ 
19 January (from Set B). Compare this to Chart IV-5 containing the correspon- 
ding analysis from Set A. The differences between these analyses are shov/n in 
Chart lV-37. Note the two regions with differences in excess of ninety (90) 
meters — on the first analysis in the sequence. Table V-6 indicates an RMS 
height difference of about fifteen meters at 250 MBS — throughout the two-day 
sequence. Chart lV-34 contains the 250 MB Height Analysis for OOOOZ 21 
January. Compare this to Cnart IV- 17 in Set A. And note the height differences 
for 21 January on Chart IV-42. Although the differences have grown in the 
eastern Pacific (to 120 meters in one place), the RMS difference has increased 
only slightly (from 13.7 to 15.S meters). 

Charts IV- 30 and lV-31 contain the (Set B) 250 MB Tempera- 
ture and Wind Analyses, respectively, for OOOOZ 19 January. Compare these to 
Charts lV-6 and IV-7 which show the corresponding baseline analyses. The (Set 
B-Set A) differences for these parameters are shown in Charts lV-40 and IV-41. 
The temperature differences are minor — but the wind differences are as great 
as fifty knots at one spot (20N 140W)! In Table V-6 one notes the RMS 
temperature difference of 0.11 degrees, and the RMS wind differences of 2.55 
and 1.78 m/sec for the u- and v-components, respectively. Thus, the wind 
differences at 250 MBS are about five (5) times greater than at 900 MBS. 

Charts IV-35 and lV-36 contain the 250 (Set B) MB Tempera- 
ture and Wind Analyses, respectively, for OOOOZ 21 January. Compare these to 
Charts IV- IS and IV- 19 which show the corresponding baseline analyses. The (Set 
B-Set A) differences for these parameters are shown in Charts lV-45 and lV-46. 
In the Pacific the temperature differences exceed one degree Celsius in three 
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regions — and two degrees in one of them. Compare IV-45 with IV-40, noting 
that the cumulative temperature differences have become increasingly signifi- 
cant. In terms of RMS, the differences doubled between 19 and 20 January. On 
Chart IV-46, the wind speed difference is 35 knots at only one location (35N 
I45W), but evident In the entire Pacific and Atlantic oceans to one degree or 
another. 
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ni. PREDICTION MODEL DEVELOPMENT 


A. Description of Models 

The basic primitive equation atmospheric prediction model employed 
in this study uses a northern hemisphere polar stereographic grid in the 
horizontal and a sigma coordinate in the vertical. This model (PECHCV) has five 
sigma layers and a 63x63 horizontal grid (381 km at 60°N). PECHCV was 
developed by ODSI in support of the SEASAT Program and was described in 
detail by Wellck (1977). A very brief description is given here. 

Conservation forms of the difference equations based on the Arakawa 
technique are integrated using a twelve-minute time step. Pressure-gradient- 
force terms are replaced by a single geopotential gradient on local pressure 
surfaces to reduce inconsistent truncation error [Kurihara (1968)]. A nonlinear 
pressure smoother is used along with momentum and temperature diffusion to 
control model behavior and computational noise. The horizontal boundary 
conditions are a persistence region below 5°N, a blend region from 5° to 20°N, 
and a fully-active region above 20°N. Centered time differencing with time 
averaging of the pressure-gradient-force term in the momentum equations is 
used. Robert (1966) time filtering of the temperature and moisture solutions is 
used for computational stability with a larger-than-usual integration time step. 

The moisture and heat source/sink terms are modeled in a manner 
similar to those in the (early) Mintz and Arakawa general circulation model as 
described by Langlois and Kwok (1969). Terms representing evaporation and 
large-scale condensation, sensible heat exchange, parameterized cumulus con- 
vection and precipitation, and solar and terrestial radiation are included. Dry 
convective adjustment precludes hydrostatic instability. Stress is applied at the 
lowest level. 

The objectives of the prediction modeling task were as follows; 
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• to implement and test an iterative procedure (based on the work 
of Tempeiton) to dynamically condition the initial state specifica- 
tion for the forecast model. 

• to modify the forecast model (initializatiwi and integration sec- 
tions) to permit data ingestion and dynamic conditioning at a 
second time origin (intermittent assimilation). 


• vary the iterative procedure to study model behavior and to 
determine a practical number of Temperton "orbits" for dynamic 
conditioning. 

• compare the behavior, efliciency, and results of each forecast 
model execution using the initial conditions being produced by the 
analysis models. 


These objectives were translated into three types of forecast model 
test demonstrations; 

• a "static initialization" baseline forecast (Chart Set E). 

• forecasts which dynamically conditioned the initial state (using a 
specified number of Temperton "orbits") prior to commencing the 
prediction, where 

Chart Set F pertains to the case in which six (6) orbits were 
made. 

Chart Sets Gl, G2 and G3 pertain to the cases in which 
twelve (12) orbits were made. 

• forecasts which dynamically conditioned the initial state (using 
twelve orbits), then integrated forward for twelve forecast hours 
and ingested new (analyzed) values at specified grid points; then 






dynamically conditioned the (new) state prior to completing the 
remainder of the forecast. These gave rise to Chart Sets HI, H2 
and H3. 

Forecasts E, F, Gl and HI were initialized with analyses for OOOOZ 19 
January. Use OOOOZ 21 January analyses for verification. Forecasts G2 and H2 
were initialized with analyses for 1200Z 19 January. Verify these at 1200Z 21 
January. Forecasts G3 and H3 were initialized using analyses for OOOOZ 20 
January. These verify at OOOOZ 22 January. 

Table V-7 contains the central processor (CP) and input/output (I/O) 
times (in seconds) for the various types of forecast runs in this study. All 
production test runs were made using a CDC CYBER 74 System. Note that the 
basic model (static initialization) takes just over an hour of CP time for a two- 
day forecast, but the I/O time is twice as large. Each Temperton orbit adds 
about 190 CP-seconds and 400 I/O seconds to these amounts. Both the model 
energetics and the dynamic conditioning of the initial state specification 
improve as the number of such orbits (taken) increases. 

1. Static Initialization 


The forecast model requires an initial state specification (of 
model variables) at each grid point in the three-dimensional lattice. Usually, 
this initial-state specification is obtained from a set of analysis model outputs or 
is derived therefrom. The analysis models output gridded values of sea-level 
pressure and sea-surface temperature, as well as the temperatures, heights and 
winds at twelve (i2) pressure levels from 1000 to 100 MBS. Heights and 
temperatures above 100 MBS are produced using regression equations. Winds 
above 100 MBS are derived from these heights. The moisture variable is derived 
from the (normalized) vorticity distribution at each pressure level using a 
procedure suggested/used by Kesel and Lewit (1974). 

Tlie forecast model requires its initial values on its terrain- 
following "sigma” surfaces [Philiips U9!»7)] rather than on pressure surfaces. 
This means that the analysis model outputs (and fields derived therefrom) must 
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be interpolated to sigma surfaces and rearranged into suitable arrays and/or 
formats. Thus, the quality of the relationships between the variables suffers 
considerably in the process. This approach to initialization shall be referred to 
as "static initialization" in this Final Report. 

Having completely specified the initial state of the atmosphere in 
terms of grid point values in the forecast model, this specification is then 
numericaiiy marched forward in time according to the primitive equations to 
obtain a forecast of the state of the atmosphere at some later time. 

2. Dynamic Initialization 

No attempt is made in the static initialization process to ensure 
that the specification of initial conditions for the forecast model results in a set 
of conditions that are numerically and dynamically consistent with the forecast 
model equations. If the mass and wind fields are not dynamically conditioned, 
gravity waves will be excited. However, when starting from a static initial 
state, it is assumed that the relationships are fairly good, and that the various 
numerical devices will bring the model trauma under control within the first six 
to twelve forecast hours with little damage to the forecast. 

An early "'pproach to specifying suitably conditlo'^-'' ,i initial fields, 
carried over from the days of filtered models, was to assume that the required 
"balaiice" could be expressed as a diagnostic relationship between the mass and 
wind fields. The simplest such relationship is the geostrophic approximation; and 
more sophisticated diagnostic relationships lead to a hierarchy of balance 
equations for the rotational component of the wind and to various forms of the 
w-equation for the divergent component. 

Although the diagnostic relationships could be used to produce a 
reasonable set of initial winds on nressure surfaces, problems cited earlier still 
had to be overcome. Specifically, we refer to: (1) the requirement for 

numerical and dynamic consistency with the forecast model; (2) the need to 
reduce initialization shock; (3) the desire to account for secondary physical 
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effects in the relationships (e.g., convection, heating, diffusion, friction). Such 
considerations have led to an alternative approach to initialization —to employ 
the forecast model itself to arrive at a set of numer». •»lly and dynamically 
consistent winds. The two main approaches to dynamic initialization are 
described below. 


In Mivakoda and Moyer's (1968) initialization scheme, the first and 
second time derivatives of the divergence are constrained to be zero, as in the 
derivation of the u) -equation; but, instead of solving explicitly a diagnostic 
equation, the primitive equations are used in an iterative fashion to alter the 
wind field slowly until the constraints on the divergence are satisfied. The 
principal computational drawback to this scheme is that it involves the solution 
of Poisson equations, repeated many times. 

In the scheme proposed by Nitta and Hovermale (1969), no 
constraints are placed on the divergence; instead, the forecast model equations 
are used to integrate backw'ards and forwards (around the initial time) using a 
scheme designed to damp the high-frequency gravity modes wl.ile retaining the 
lower-frequency balanced meteorological modesi. Besides lifting the restrictions 
imposed by the Miyakoda-Moyer scheme, the Nitta-Hovermale scheme is more 
adaptable since provision can be made for the mass field to adjust to the wind 
field. The main disadvantage of the scheme is the somewhat slow rate of 
convergence towards the required solution. 


The dynamic initialization scheme used in this study is a variation 
of the Nitta-Hovermale scheme proposed by Temperton (1973), and implemented 
and tested in a multi-layer (PE) model by Kesel and Wellck (1975). One iteration 
(or "orbit") of this scheme is carried out as follows; beginning at a central time 

level, T , integrate forward to (x + NAt), and then integrate backward from t 

0 0 o 

to ( r - NAt). To obtain the new estimate of the wind fields at time level x , 
restore the mass fields and replace the wind fields by 


and 


u* = [vi( i^^+NAt.) - ii(l^^-NAt;)] 
v' = 'a [v(v^+NAt) + v(x^-NAt)3 
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This process may be repeated for "N" orbits until successive wind sets differ 
from one another by less than epsilon. As a practical matter, we have chosen to 
use the layer-mean squared mass divergence as the measure for assessing the 
quality of the conditioning. For this, we define mass divergence as 

M ^ (4>u) + Ip- (<lv) 

where < 1 ? is the geopotential-height. 

A version of the northern hemisphere primitive equation atmos- 
pheric forecast model PECHCV which includes Temperton-style dynamic initiali- 
zation was constructed for this study. This version of the forecast model, known 
as DYNMOD, has provision for a varying number of orbits of variable length to 
be performed (after the static initialization) in order to dynamically condition 
the initial state. The orbit length used for this study was NAt = 1 hour or five 
time steps. Thus, each orbit was computationally equal to two forecast hours 
(roughly). 


3. Data Assimilation Method; Dynamic Conditioning at Two Time 

Origins 

A data assimilation experiment using the DYNMOD forecast 
model was designed to simulate the effect of ingesting off-time satellite data on 
a forecast. This experiment combines dynamic initialization with the (inter- 
mittent) assimilation of data during the course of the forecast. The steps 
involved in the data assimilation experiment are listed below: 

• static initialization 

• dynamic balancing with twelve orbits 

• the first twelve hours of the forecast 

« ingestion of mass field data over a limited geographic region 

• dynamic balancing with twelve orbits 

• continuation of the forecast to ^8 hours 

Data was ingested into the forecast model in a pie-shaped wedge 
over the Pacific extending from 130W to 155E. In order to avoid boundary 
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discontinuities created by t'.ie data ingestion process, the 75° wedge was divided 
into two blend zones of 25° each and a center zone in which data was directly 
substituted. Only mass field data (sea level pressure, temperatures and 
geopotentials) were ingested. The winds from their twelve-hour forecast values 
were forced to adjust during the dynamic conditioning that followed the data 
ingestion. 


B. Test Results 


Three types of forecasts were made in the test demonstration task; 
(1) static initialization; (2) dynamic initialization using a varying number of 
Temperton orbits; (3) dynamic initialization at two time origins with data 
assimilation. Selected charts will be introduced and dis-^ussed in order to allow 
the reader to determine the validity of the forecasts, as well as the comparative 
performance. Graphs that contain plots of sensitive model parameters (RMS 
mass divergence; kinetic energy; number of precipitating grid points) for 
forecasts initialized with either zero, six or twelve Temperton orbits will be 
presented. 


1. Static Initialization (Forecast E) 

Forecast Run E (Static Initialization) was initialized using analysis 
model outputs for OOOOZ 19 January. Charts IV-I and IV-5 contain the starting 
sea-level pressure and 50f MB height analyses, respectively. The verification 
analyses for OOOOZ 21 January are provided as Charts IV-12 and IV-16, respec- 
tively, for the SLP and Z500 analyses. 

The Run E outputs are provided as Charts IV-73 and IV-74, the 48- 
hour SUP and Z500 forecasts, respectively. Run E deepened the low near Japan 
from 976 MBS to 960 MBS, and moved it to 53N 180. It actually verified at 956 
MBS slightly north of the forecast position. Run E also deepened the low near 
Newfoundland from 996 MBS to 988 MBS and moved it northeastward to 44N50W. 
It verified about 350 miles further north and slightly deeper than predicted (at 
984 MBS). An intense low south of Novaya Zemyla (984 MBS) was predicted to 
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deepen slightly to 980 MBS, but it actually filled a little to 988 MBS. Over the 
United States, Run E predicted a 1004 MB low over Illinois — but it verified in 
northern Mississippi at 992 MBS. 

Error statistics for Run E were encouraging, however. The RMSE 
was 4.99 MBS in the sea-level pressure forecast. Compare this to an RMS actual 
change (between the initial and verification analyses) of 8.36 MBS during the 48- 
hour period. At 500 MBS, Run E had an RMSE of 46.4 meters compared to an 
RMS actual change of 69.7 meters. Thus, Run E exhibited good skill when 
compared to persistence. 

2. Dynamic Initialization (Forecasts F and Gl) 
a. Geostrophic Adjustment Study 

The prediction model DYNMOD (which stands for "dynamic 
initialization model") was used to produce 48-hour forecasts after first complet- 
ing six (6) Temperton orbits (Forecast Run F) and twelve (12) Temperton orbits 
(Forecast Run Gl). Each of these forecasts may be compared to the Run E 
forecast discussed earlier, since all three used the same initial conditions (for 
OOOOZ 19 January ). 


To measure the degree of geostrophic adjustment (dynamic 
conditioning) in each forecast, statistical measures of sensitive model param- 
eters were produced each forecast hour. Perhaps the most useful measure is the 
RMS mass divergence at each model level. Figures VI- 1 through VI-5 show the 
effect of "orbital integrations" at each of the five model levels, respectively. 
On each figure, there are three curves; Curve A depicts the static initialization 
forecast (Run E); Curve B depicts Run F which employed six Temperton orbits. 
Curve C depicts Run Gl which employed twelve orbits (before commencing the 
forecast). 


From an examination of these figures, it is obvious that 
dynamic initialization has a beneficial effect on the balance between the mass 


and motion fields estimates. (It remains to be seen if the forecasts are any 
better.) At the sigma=0.9 level (shown in Figure VI-1), the mass divergence 
parameter increases from five to twelve units (on Curve A) within the first 
twelve forecast hours — but settles back to about six units toward the end of the 
forecast. Curve B (six Temperton orbits) shows less trauma. And Curve C 
(twelve orbits) shows even less trauma. Observe, also, that the "envelope” 
described by the three curves "narrows" appreciably after the first 24-36 
forecast hours. This particular effect (narrowing of the envelope) is even more 
pronounced at other model levels. The suggestion is that the model itself has the 
adjustive power to accomplish the job, but it would take the model longer! One 
detects a slight computational "edge" to the task of accomplishing ein acceptable 
degree of geostropiiic adjustment through the iterative process being tested, 
rather than let the model do it at the same time the forecast solution is 
evolving. 


Thus, our experience compels us to say that dynamic condi- 
tioning of an initial state specification can be substantially accomplished before 
the forecast commences, thereby eliminating or reducing harmful, spurious 
physical effects. 


Recall, however, that two general types of Temperton orbits 
can be employed. The first type has been demonstrated herein. The second type 
has not been discussed. In the second type, the winds are restored and the mass- 
field variable estimates are averaged at the end of each orbit. From geostrophic 
theory, both types probably should be done — perhaps in an alternating manner. 
And the two solutions could be combined according to the latitude (since the 
winds oscillate about the mass-field gradients poleward of about 30° latitude, 
and vice versa equatorward of 30° latitude. The second type of orbits, however, 
was not tried in this study. 

Figures VI-6 through VI- 10 show the effect of initialization 
(type) on the layer-mean kinetic energy. These show that dynamic initialization 
has far less effect on mean kinetic energy (than on mass divergence), which is as 
expected. The "envelopes" tend to be quite narrow in these figures as well. 
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Even so, by taking the orbits the model commences the forecast (generally) with 
KE values closer to the values the forecast ultimately produces. 

Figures VI-1 1 ^d VI-12 show the effect of initialization on 
the precipitation mechanisms in this model. Two types of precipitation are 
modeled; large-scale and convective. Note that here (too) the envelope narrows 
as the forecast lengthens. Thus, the impact on convective precipitation is to 
reduce the affected areas and amounts in the first forecast day, with little 
impact after the adjustment has matured. In Figure VI- 12, however, the effect 
of dynamic initialization is to increase the large-scale precipitation area and 
amounts by about 10-20%. A reasonable explanation is that the vertical velocity 
patterns become rather well developed in the orbital integration period. Im- 
provements in modeling of "secondary" physical effects such as precipitation 
may be one of the more important motivations to dynamically condition the 
initial state of a (PE) prediction model. 

b. Chart Discussion 
(1) Forecast Run F 

Run F was dynamically initialized with only six Temper- 
ton orbits. Even so, the alterations in the wind fields were large. The resulting 
forecasts should be compared to Run E (static initialization) and/or to Run G1 
(twelve orbits). 

Charts IV-75 and IV-76 show the 48-hour se. • I'-vel pres- 
sure and 500 MB height forecasts (for Run F), respectively. The corr; spending 
verification analyses are provided a;; Charts IV-12 and lV-16. 

Run F tended to produce lows which were 1-7 millibars 
weaker than in Run E. The 3apan low was 1-2 millibars weaker than in Run E. 
With the Asian low Run F was actually better than Run E by not deepening it as 
much. The differences between the two forecasts for both the U.S. and 
Newfoundland lows were negligible. At 500 MBS, Run F and Run E look almost 
identical. If anything, the major low centers appear to be slightly weaker in Run 
F. 
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Statistically, Run F was slightly better than Run E. At 
sea-level, the Run F RMSE was 4.97 millibars compared to 4.99 for Run E. At 
500 MBS, Run F was 46.2 meters compared to 46.4 meters. Strictly speaking, 
the differences are not statistically significant. 

(2) Forecast Run G 1 


R in G1 was dynamically initialized using twelve (12) 
Temperton orbits. When the results are compared to other forecasts (E,F) in the 
same initialization class, it is difficult to see any significant differences. The 
pressures of the major storm centers are identical to Run F. And, as with Run F, 
Run G1 seems to predict lows which are slightly weaker than in the static 
initialization forecast (Run E). The same generalization holds at 500 MBS, where 
it is difficult to detect differences between Runs E, F and Gl except to say that 
the centers of lows seem to get progressively weaker as the number of orbits is 
increased (before commencing the forecast). 

Charts lV-77 and lV-78 contain the 48-hour sea-level 
pressure and 500 MB height forecasts, respectively. Compare these to Charts 
IV-75 and lV-76, respectively, for Run F; and to Charts IV-73 and lV-74, 
respectively, for Run E. 


Statistically, Run Gl produced RMSE values nearly iden- 
tical to those produced by Runs E and F. 

In summary, it has been demonstrated that the quality of 
the mass-motion parameter relationships is improved as the number of orbits is 
increased. Further, one detects an impact on the large-scale precipitation 
mechanism. Yet, the differences between the static and dynamic initialization 
two-day forecasts were not significant. 
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3. Data Assimilation; Dynamic Conditioning at Two Time Origins 


a. Run HI vs. Run Gl 


In Forecast Run HI, the sequence of events was as follows: 
(1) static initialization; (2) twelve orbits; (3) forward integrations for twelve 
forecast hours; (4) data assimilation in a 75° (longitude) band in the central 
Pacific; (5) twelve orbits; (6) forward integrations for last thirty-six hours. 

Run HI was initialized with OOOOZ 19 January data, a time 
at which the data base is best (generally) in the Pacific area. The ingested data 
came from an analysis for 1200Z 19 January -- a time at which the data base in 
the Pacific is not very rich. Thus, we may have introduced a problem by 
updating our forecast solution with values that may not be as valid as those in 
the forecast solution. If this is so. Run H2 may be helpful -- since it was 
initialized at 19/1200Z and was updated using 20/0000Z analysis values. 

The 48-hour sea-level pressure and 500 MB height forecasts 
for Run HI are shown in Charts IV-79 and IV-SO, respectively. The verification 
analyses are provided as Charts IV- 12 and IV'- 16. They may also be compared to 
Forecast Run Gl which was dynamically initialized with twelve orbits also — but 
which did not assimilate new values at a second time origin. The Gl charts are 
lV-77 and lV-78. 


The biggest difference between Gl and HI is in the western 
Pacific — where HI produced a low of about 966 MBS, whereas Gl deepened it to 
961 MBS. For the other major features. Runs HI and Gl produced nearly 
identical solutions. The RMS statistics support this assessment. The 48-hour 500 
MB height RMSE values differed by only 0.1 meter; and the sea-level pressure 
RMSE values differed by only 0.03 millibars. Thus, the only major difference 
occurred at the upstream edge of the assimilation region. 
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b. Run H2 vs. Run G2 


Another pair of forecasts was initialized using 19/1200Z 
analyses. These are Runs G2 and H2. G2 employed twelve orbits at the initial 
time origin only. H2 was also conditioned this way at the initial time origin, but 
assimilated updated values in the 75° band in the central Pacific twelve hours 
into the 48-hour forecast. In this comparison, H2 did produce a slightly better 
forecast. At 500 MBS, H2 had an RMSE of 49.2 meters while G2 had 50.1 
meters. At sea level, H2 had an RMSE of 5.21 millibars while G2 had 5.26 
millibars. 


The SLP and Z500 48-hour forecasts for Run H2 are provided 
as Charts lV-83 and IV-84, respectively. Compare these outputs to Run G2 
outputs contained in Charts lV-81 and lV-82. (The initial conditions for both runs 
are shown in Charts lV-8 and IV-9. The verification analyses are provided as 
Charts lV-23 and lV-24.) 

Neither H2 nor G2 handled (well) the Japan low (which moved 
to the Aleutians during the period). H2 filled it from 972 MBS to 974 MRS, even 
though it deepened to 960 MBS by 21/1200Z. Run G2, in contrast, kept the 
central pressure about the same. Neither forecast deepened the system. 

During this period, a major low evolved in the eastern U.S. 
Roth H2 and G2 produced 998 MB lows near Chicago. The system actually 
deepened to 984 MBS.... from a 1007 MB low near northeastern Nebraska on 
19/1200Z. 


The 980 MB low (at 19/1200Z) just south of Newfoundland 
filled to about 1000 MBS during the period. Both G2 and H2 only filled it to 992 
MRS (about half as much as it did in fact). 


Although the H2 forecast was slightly better than the G2 
forecast (statistically), one is hard pressed to detect it with an examination of 
the output charts. 
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c. Run H3 vs. G3 


Here are two more runs from the "G" and "H" series. (Recall 
that data assimilation takes place in the "H" series.) Runs H3 and G3 were 
initialized using analyses for OOOOZ 20 January (Charts IV- 10 and IV- 11). The 
verification analyses are provided as Charts IV-25 and IV-26. 

Once again, H3 appears to be better statistically. Its RMSE 
at 500 MBS was 48.6 meters, compared to 50.0 meters for Run G3. At sea level, 
H3 had 5.27 MBS compared to 5.33 MBS for G3. 

The H3 forecasts are shown in Charts lV-87 and IV-SS. 
Compare these to the G3 outputs in Charts IV-85 and lV-86. 

Clearly, the only difference at sea level takes place in the 
northwest Pacific, where H3 produces a 970 MB low and G3 produces a 974 MB 
low. It actually verified at 968 MBS. 

In the eastern U.S., both H3 and G3 predicted the low to be 
988 MBS — and it verified about ten millibars deeper than that. 

4. Comparative Performance Summary 

Table I1I-8 shows the error statistics for all production test runs. 
The most striking aspect of the statistics is that the differences are quite small. 
For all four forecasts made from 19/OOOOZ, for example, the RMSE ranges were 
from 4.97-5.01 millibars at sea level, and from 46.2-46.4 meters at 500 MBS. 
For the two forecasts initialized at 19/1200Z, as well as for the two forecasts 
initialized at 20/0000Z, one can detect slight superiority in the H2 and H3 
forecasts. 


Table 111-9 shows comparative model performance at sea level for 
the four forecast runs initialized at 19/OOOOZ. The central pressures for four 
systems are tabulated. The entries in parentheses represent our best estimates 
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of the actual central pressures, whereas the associated entries (without paren- 
theses) indicate the value of the last closed isobar. The 3apan low (which moved 
to the Aleutians) is most relevant to this study -- since data ingestion occurred 
there (in Run HI). For Runs F and Gl, which were initialized with six and twelve 
orbits, respectively, the low appears to be slightly weaker (1-2 MBS) than in Run 
E (static initialization). In Run HI, dynamic conditioning at two time origins 
(with data ingestion, too) produced a low which was even weaker. With the other 
features, the differences were extremely minor. But, we still detect a slight 
"flattening" of the pressure profiles through the major lows when dynamic 
initialization takes place. 

On the basis of the small number of forecasts made, we are 
hesitant to draw strong conclusions about the benefits of dynamic initialization 
and/or data assimilation. On the basis of these tests it is tempting to conclude 
that the impact of dynamic initialization on the forecasts is minor — even 
though it improves the quality of the coupling in the initial state specification 
and impacts on secondary physical effects (the number of precipitation grid 
points, for example). 


Data assimilation was tried in three forecasts (HI, H2, H3). We 
detect slightly better statistics at 500 MBS. In the sea level pressure forecasts 
we are not sure we detect any improvement in the forecasts. Since the 
computational burden was so much greater in the G-series and H-series fore- 
casts, one must carefully consider the cost-to-benefit ratio. Perhaps the data 
assimilation issue can (only) be resolved when large quantities of accurate, 
representative data are available for the assimilation. It doesn't make sense to 
ingest values from a later analysis — if such an analysis is based on an equally- 
poor data base. 
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CHART IV-3: 900MB Wind Analysis, OOOOZ 19 January 1979. ^rt Set A, 












CHAF<.T IV-6: 250MB Temperature Analysis, OOOOZ 19 Oanuary 1979. Chart Set A. 
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CHART 500MB Height Analysis, J200Z 19 January 1979. Chart Set A. 












CHART IV-12; Surface Pressure Analysis, OOOOZ 2! lariuary 1979. Chart Set A 
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CHART lV-13: 900MB Height Analysis, OOOOZ 21 3aiuary 1979. Chart Set A. 






CHART IV- 14: 900MB Temperature Analysis, OOOOZ 21 January 
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CHART IV-27: 900MB Temperature Analysis, OOOOZ 19 January 1979. Chart Set B. 
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CHART 1V-3I: 250MB Wind Analysis, OOOOZ 19 Oanuary 1979. Chart Set B, 











CHART IV-33: 900MB Wind Analysis, 00002 21 January 1979. Chart Set B. 





CHART lV-35: 250MB Temperature Analysis, OOOOZ 21 January 1979. Chart Set B. 
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CHART IV- 36: 250MB Wind Analysis, OOOOZ 21 January 1979. Chart Set B. 
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CHART IV-39: 900MB Wind Difference, OOOOZ 19 January 1979. Chart Set B minus Chart Set A. 




CHART IV-41: 250MB Wind Difference, OOOOZ 19 January 1979. Chart Set B minus Chart Set A. 



CHART lV-42; 250 MB Height (Cumulative) Difference, OOOOZ 21 January 1979. Chart Set B minus Chart Set A 







CHART IV-46: 250MB Wind (Cumulative) Difference, OOQOZ 21 January 1979. Chart Set B minus Chart Set A 
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CHART IV-48: 900MB Temperature Analysis, OOOOZ 19 January 1979. Chart Set C. 
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CHART lV-49; 900MB Wind Analysis, OOOOZ 19 3anuary 



CHART IV-50; Surface Pressure Analysis, OOOOZ 21 January 1979. Chart Set C 





CHART IV-51: 900MB Temperature Analysis, OOOOZ 21 January 1979. Chart Set C 





CHART JV-52: 900MB Wind Analysis, OOOOZ 21 January 1979. Chart Set C 






CHART lV-53: Surface Pressure Difference, OOOOZ 19 January 1979. Chart Set C minus Chart Set A. 
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CHART IV-59: Surface Pressure Analysis, OOOOZ 19 January 1979. Chart Set D. 
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CHART IV-60: 900MB Temperature Analysis, OOOOZ 19 January 1979. Chart Set D, 


















CHART IV-63: 900MB Temperature Analysis, OOOOZ 21 January 1979. Chart Set D, 
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CHART IV-64: 900MB Wind Analysis, OOOOZ 21 January 1979. Chart Set D. 



CHART lV-65: Surface Pressure Difference, OOOOZ 19 January 1979. Chart Set D minus Chart Set A. 
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CHART IV-68: Surface Pressure (CumuJative) Difference, OOOQZ 21 January 1979. Chart Set D 
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CHART IV-70: 900MB Wind (Cumulative) Difference, OOOOZ 21 January 1979. Chart Set D minus Chart Set A, 





CHART IV-71: Surface Pressure Difference, OOOOZ 19 January 1979. Chart Set D minus Chart Set C. 







(Cumulative) Difference, OOOOZ 21 January 1979. Chart Set D minus Chart Set C 



CHART lV-73; 48-Hour Surface Pressure Forecast from OOOOZ 19 January 1979. Chart Set E. 











CHART lV-75; 48-Hour Surface, Pressure Forecast from OOOOZ 19 January 1979. Chart Set F, 




CHART IV-76: 48-Hour 500MB Height Forecast from OOOOZ 19 January 1979. Chart Set F, 



CHART IV-77: 48-Hour Surface Pressure Forecast from OOOOZ 19 January 1979. Chart Set G1 









CHART IV-80: *»8-Hour 500MB Height Forecast from OOOOZ 19 January 1979. Chart Set HI. 
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CHART fV-82: 48-Hour 500MB Height Forecast from I200Z 19 nan'!?'’y 1979. Chart Set G2, 

















CHART IV-85: 48-Hour Surface Pressure Forecast from QOOOZ 20 January 1979. Chart Set G3 






CHART IV-86: ifS-Hour 500MB Height Forecast from OOOOZ 20 January 1979. Chart Set G3. 











TABLE V-I: RMS DIFFERENCES BETWEEN "SET A" AND "SET C" ANALYSES 


PARAMETER 

DATE* 

3AN 19 

OAN 20 

3AN 21 

PS (mbs) 

0.62 

0.67 

0.55 

AZ900 (m) 

2.06 

3.43 

3.00 

AT900 (°C) 

0.02 

0.04 

0.03 

ALI900 (m/sec) 

0.46 

0.71 

0.55 

AV900 (m/sec) 

0.55 

0.66 

0.53 

AZ250 (rn) 

0.96 

2.37 

1.92 

AT250 (°C) 

0.02 

0.03 

0.05 

AU250 (m/sec) 

0.18 

0.39 

0.27 

AV250 (m/sec) 

0.22 

0.40 

0.28 


*At 0000. 
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TABLE V-2; RMS HIFFERENCES BETWEEN "SET A" AND "SET D" ANALYSES 


PARAMETER 

DATE* 

JAN 19 

JAN 20 

JAN 21 

1 

PS (mbs) 

O.Sl 

O.Sl 

0.76 

AZ900 (m) 

3.23 

2.7S 

2. 84 

AT900 (°C) 

0.03 

0.04 

0.04 

AU900 (m/sec) 

0.63 

0.64 

0.60 

AV900 (m/sec) 

0.76 

0.65 

0.63 

AZ250 (m) 

1 . 64 

1.66 

1.76 

AT2-50 ("O 

0.02 

0.03 

0.03 

AU250 (m/sec) 

0.21 

0.2S 

0.28 

AV250 (m/sec) 

0.32 

0.30 

0.27 


♦At 0000 Z. 
































TABLE V-6; RMS niFFERENCES BETWEEN ’’SET A” AND ’’SET B” ANALYSES 


PARAMETER 

DATE* 

JAN 19 

JAN 20 

JAN 21 

PS (mbs) 

0.00 

0.11 

0.10 

^Z900 (m) 

l.OS 

1.31 

1.16 

AT900 (°C) 

0. 17 

0.32 

0.29 

AU900 (m/sec) 

0.65 

0.63 

0.7S 

AV900 (m/sec) 

0.30 

0.32 

0.31 

AZ250 (m) 

13.7 

17.5 

1.5.S 

AT250 (‘-’’O 

0. 1 1 

0.22 

0.26 

AU250 (m/sec) 

2.55 

2.97 

1.70 

AV250 (m/soc) 

1.7S 

2.72 

1.90 









TABLE V-7: COMPUTER EXECUTION TIMES FOR TWO-DAY FORECASTS 
USING THE CDC CYBER 74 SYSTEM. 


Model/Run Type 

Computer Time (seconds) 

CP 

I/O 

• Run E 

Static Initialization plus 2-day forecast 

3,770 

7,460 

• Run F 

Six Orbits plus 2-day forecast 

4,900 

9,874 

• Runs GI, G2, G3 

Twelve Orbits plus 2-day forecast 

6,020 

12,200 

• RunsHl,H2, H3 

Twenty-four orbits w/ data assimilation 

8,340 

17,000 


Note; Each Temperton "orbit" takes about 190 CP-seconds and about 
400 seconds of I/O. 









TABLE V-8; COMPARATIVE ERROR STATISTICS FOR 48-HOUR FORECASTS 


Starting Forecast Sea- Level 

Time Identifier Mean 



G2 0.02 

19/12Z 

H2 0.08 


G3 -0.01 

20/00 Z 

H3 0.04 


19/OOZ H4 -0.17 


Pressure I 300 MB Height 


RMSE 

Mean 

RMSE 

4.99 

1 

-2.62 

46.4 

4.97 

-2.52 

46.2 

4.98 

-2.04 

46.3 

5.01 

-1.08 

46.4 

5.26 

-1.19 

50.1 

5.21 

-0.65 

49.2 

5.33 

3.55 

50.0 

5.27 

3.U2 

48.6 


4.99 -2.04 46.2 























TABLE V-9; COMPARATIVE FORECAST MODEL PERFORMANCE FOR 
FOUR TYPES OF INITIALIZATION 



*48-Hour Forecasts 



VI. FIGURES 


VI- 1 







:NCE at SI6MA=0.7 level 






FIGURE VI-4; EFFECT OF INITIALIZATION ON RMS MASS DIVERGENCE AT SIGMA=0.3 LEVEL 



FIGURE VI-5: EFFECT OF rilTIALIZATIOfi ON RMS MASS DIVERGENCE AT SIGMA=0.1 LEVEL 



FIGURE VI-6: EFFECT OF INITIALIZATION ON LAYER MEAN KINETIC ENERGY AT SIGMA-0.9 LEVEL 




FIRUPE '71-7; EFFECT OF INITIALIZATION ON LAYEP MEAN KINETIC ENERGY AT SIGMA-0. 7 LEvEL 




FIGURE VI-8; EFFECT OF INITIALIZATION ON LAYER MEAN KINETIC ENERGY AT SIGMA=0.5 LEVEL 





FIGURE VI-9; EFFECT OF INITIALIZATIOM ON LAYER MEAN KINETIC ENERGY AT SIGMA=0.3 LEVEL 




FIGURE VI-10; EFFECT OF INITIALIZATION ON LAYER f€AM KINETIC ENERGY AT SI6MA=0.1 LEVEL 





GDRf VI-12; EFFECT OF rilTIALlZATION 0*: 1 ARG^: GCALF PkECIP I TATiOf! MECHANIGM 
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APPENDIX A 


SCALAR ANALYSIS USING THE PATTERN -CON SERVING TECHNIQUE 
I . INTRODUCTION 

The goal of the PCT analysis is to blend the following 
information: the new data; the most recent past analysis 

(or forecast) value (the first guess); the gradients (of 
the first guess) in eight directions from each grid point; 
and the Laplacian (of the first guess). The relative impor 
tance of each piece of information is specified by an array 
of weights. 

The desired blend is realized by minimizing the sum 
of the deviations of the various characteristics of the 
analysis from their counterparts in the first guess. The 
minimization is accomplished with an elementary application 
of the calculus of variations. 

Information is spread spatially by the gradient and 
Laplacian terms. In a surface analysis, there are some- 
times natural obstacles (mountain ridges, coastlines, etc.) 
beyond which an analyst would not allow a new observation 
to influence the analysis. This kind of constraint can 
be simulated in the objective analysis by reducing the 
weights of the gradients and Lap'lacian along the demarca- 
tion zone. 
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An analysis cycle consists of three major steps: 

• Assemble the data at grid points. 

• Solve the mi n imi zat ioi. equation. 

• Re-evaluate the weight of each report. 

In order to adequately evaluate the weight of each report, 
at least two cycles are required. 


II. ASSEMRLY 


We shall refer to the guess field as ^ with weight 
A. .. On the first cycle, it is the first guess, and A. . 
has a low and probably uniform value. On subsequent cycles, 

P. . is the result of the previous cycle, but A. . keeps 

» j ^ * j 

its original value. 

The purpose of the assembly procedure is to incorporate 
the observational data into the first guess field P. ., 
taking into account the subjective specification of each 
report's reliability (DWT) and its distance from the grid 
point. Grid points within a specified influence region of 
each observation are affected by that observation. The size 
and shape of the influence function are determined by the 
data density, analysis cycle number and first guess field shape 
(i.c., gradient and Laplacian) , respectively. An information 
density field is used to produce a factor (FACT) which varies 
the basic radius of the influence for each observation between 
a minimum and maximum limit. In areas of dense data concen- 
tration, the influence radius is set to the minimum value 
so as not to spread a data report's influence so far that it 
interferes with the already well-specified observed values. 
However, if the observation is isolated, its influence is 
spread to the maximum allowed. 
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The assembly radius from an observation which includes 
all gridpoints to be influenced by that observation is cal- 
culated as; 


RADIUS 


FACT * AMAP * RAD 
AMESH 


where AMAP = map factor 


AMESH = standard meshlength of the reference 
latitude 

RAD 5= a multiple of AMESH 

FACT = factor proportional to the information 
density 


FACT is computed as follows: 

FACT = RADMAX - INFOFAC * (RADMAX - RADMIN) 

(I, J) 

where RADMAX = maximum allowable factor 

RADMIN = minimum allov/able factor 

INFOFAC ( I, J) = value of information density 

factor nearest observation location 

The maximum allowed radius (RADMAX) is decreased with each 
cycle in order to better define progressively smaller 
scales.... in the manner used in Cressman analyses. 
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The basic influence function has a weight of on® at 
its center (observation location), decreasing to zero at 
its maximutn radius as determined by the information den- 
sity. The fraction of the radius to which the weight value 
remains one (FRAD varies between a minimum and maximum 
value determined by the curvature or gradients of the field. 
In systems such as cyclones, the curvature and gradients 
are large and an observation’s full influence should not 
extend far from its location since it is less representative 
of the rapidly varying field. In anticyclones, these 
characteristics vary less rapidly, and it is acceptable 
to have the full weight of the observation included in the 
assembled fields at larger distances. 


FRAC is computed as follows: 


FRAC = 1 

.0 - SF 

^ , GRAD . 

where SF ^gRADMX^ 

GRAD 

= maximum gradient at 

GRADMAX 

= maximum gradient fo 

WTLAPL 

= iLaplacian ( I , J) | 

WTLAPT.M 

= Percentage of the m 
the entire field. 


' WTLAPL £ WTLAPLV 
GRAD < GRADMX 
FRACMTN < FRAC < FRAC^AX 


. WTLAPL . 
'WTLAPLM ' 


whichever 
is greater 


a gridooint 
r the entire field. 


aximum Laplacian for 
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As tho first step in the assembly procedure, the guess 
field is interpolated at the observation location and the 
difference between the observation and the guess field 
determined (DIF) . If DIF is greater than the gross tolerance 
for tho parameter being analyzed, it is excluded from the 
assembly process in that analysis cycle. It may be included 
in some subsequent cycle (s). 

Next, the value of the influence function appropriate to 
the distance of the grid point from the observation is com- 
puted (w) , where: 

- _ . distance < 

" ' -radTus- 

, - distance 
otherwise w = ^ ~ RADIUS 

1.0 - FRAC 

(ror the upper air analyses, FRAC is set to a constant.' 


For each gridpoint affected by the observation, a cumu- 
lative sum of the product (W*DWT) *W*DIF is computed at the 
appropriate I,J. Also, a field of the product W*DWT is 
accumulated. Once all observations have been processed, the 
assembled value is obtained by dividing the two fields at 
all grid points: 


NOBS 

i: [W*DWT(K) ] *W*DIF 

p = p t i^;i 

NODS 

Z W*DWT(K) 

K=1 


for 1=1, M 
J=1,M 


= assembled field value 
= number of observations 

= data weight assigned to an observation 


A-h 


P 

NODS 

nWT 


Ill 


MINIMIZING THF ERROR FUNCTIONAL 


TARLE A-l; PCT SCALAR CONSTRAINTS 



Constraint 



Weight 

j 

= Variable being analyzed (assembled 

value) 


^i, j 


= y axis gradient = P. ... - P. . 

1/1+1 i/D 

(computed from non-assemblt d value 

of first 

guess) 

®i, j 

j 

= X axis gradient = P- - P. . 
(computed from non-assembled value 

of first 

guess) 

^i, j 

j 

= x-l,yfl gradient = " Pi,j 

(computed from non-assembled value 

of first 

guess) 

^i, j 


= x+l,y+l gradient = - Pi,j 

(computed from non-assembled value 

of first 

guess) 

F . . 

1,1 

j 

= Laplacian - + P^^^. 

(computed from non-assembled value of first 

. - 4P. . 

-1 1,1 

guess) 

D . . 

1,1 


The first guess shapes y, v, a, 0 and L and their re- 
spective weights B, C, E, F and D have a constant value 
during the entire analysis . Within limits specified by the 
weights, we require the final analysis to have similar values 
of the constraints as the first guess field. 
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To effect this matching, we shall minimize the following 


integral : 


I = // r 


(a) 

'‘i.j “T.j - ‘■i.i’’ *■ 


(b) 

“i,j ‘n.j.i - ‘■i.j - * 


(c) 

5. . , (P* . _ p* . ^ _ u 

■ 4 - 

(d) 

^ (P*^i • -P^ • - V. + 

i*j i+l,] i»j i»3* 


(e) 

"’i.j - 

■f 

(f) 

* 


Cg) 

’^i.j ■ ^i4.i,j-i ' “i+1, 


(h) 

'’i.j “-Ui.:*! - n.i - - 


(i) 

"-2.3 - 

* 

(j) 

“i.j <'‘iii,j * '■t-1,3 * M.jn ^ 

P* . , - 
i»3-l 


A- 8 


Jdxdy 




In the above, the starred quantities are the analysis 
values we are seeking. Each term is a departure from the 
desired matching of differential properties. Extra terms 
have been added to account for the effect of a changing 

j on the differential properties computed at surrounding 
points. Their effect is to more closely couple neighboring 
grid points. See Figure A-1 for a depiction of the minimi- 
zation stencil as it relates to the terms of equation [A.l]. 
To minimize the integral, we simply take the first variation 

with respect to ., and set it to zero (see equation [A. 2]) 

t 3 

The solution of the resulting equation will be the PJ ^ that 
will cause the integral to be minimized. The fact that each 
term is squared ensures a minimum as opposed to a maximum 
value . 
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© 


© 


© 



LEGEND; ( ) = constraint from Equation [A.l] 

o = difference 

> = gradient 

/= laplacian at grid point 
= grid points 


FIGURE A-1: 


SCALAR MINIMIZATION STENCIL 
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- // t (Pl.j - 


J 


-2 Bi,j in.i.i - n,2 


"i.i’ 


^ 2 (P*^ . - ?l-_, - u, ,.,) 

- 2 - v,^.) 

" 2 (P*^. - 


- 2 - P-^ . 


“i,j’ 


IA.2] 


Oft* 


* 2 (P*^j - Pi+lJ_l - =i + l,j-l> 

- 2 Pi,j (PI.l,3n - n.j - 

* 2 (P‘,, - n-in-i - Si-l.j-l’ 


%y. 


4 Pt . - I-i ,) 


] dxdy 0 
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The terms in can be grouped into three categories 

1. Those involving .. 

i/D 

2. Those involving P* at surrounding points. 

3. Those not involving P*. 


S. .P* 
1/3 i 


r ‘ *i.j ’ °i,i * =i.j-i * =i.j 
■’ I * * '■i.j * '■i-l.i-l * 16 "i.j I 


A- 11 


-H 


®i,j + l " ^i.j 


^i- 1,3 '* 1 - 1 , j 


i»D 


^1-1, j+1 " ^i+1,3-1 “ ^i,j + i 

^i-l,j-l *’i-l,j-l ‘ ■* °i,j **1+1,3 ■ *’ ^i-lO 

^ °i,3 n.:.i • ^ °i,3 **1.3-1 


■G. . 
i» D 


-A. . P. . + B. . U- . - B. . , u. , + C. .V. 

J-»3 1,3 1.3 1.3 1.3-1 1.3-1 1.3 1.3 


■ ‘^i-l.j '*i-l.j 


^i,3 °i.j ” ^i+l,j-l “i+l,j-l 


+ F. . S. . - F. , . , S. , . , + 4 D. . L. . 
1.3 1.3 1-1, 3-1 1-1, 3-1 1.3 i.3 


Note that all terras in S and G except A. . in S. . and 

1 . j ^ . J 

-A. . P. . in G. . involve first-guess pattern information 

1.1 1.3 1.3 

which is consistent during the analysis. 

The minimization may be written as 


1.3 


P* 

1 


- (G, 


1.3 


H. . ) = 
1.3 


(A. 3] 


In H. . , 
1.3 

at each point, 


let us group together the coefficients of P* 






-H 


Uj 


* ‘-^l-lj - ^ Di.j> n-l.i 

* ^ °i.j» '’1^1,1 

* n-x.^.i 

. (-B, . . 4 

- (-Vi.j-r n-i.i-i 

* <-«i,j-i - -« ''i.j’ Pi.j-i 


Define 


^i.3 ' 


Y, : B. . 

i.J 1.3 

::. . : -F, 

1.3 1.3 


*^i,3 ■ .3 


Note that X, Y, Z and R have a constant va?"e during 
the analysis. 

Then : 


-H = x p* -y 

1.3 1-1,3 1-1. 3 i.3 ‘i'*-1.3 

>R. P* -Y P* 

‘ 1-1,] 1 "'i, j 1,3+'. 

^ ’i,3 '’i+Uj + l ^ ^i-L,j-I ‘’'i-l,j-l 

^i,j-l ^i,3-l '’i + l,3-l 

- ^ °i,3'^t-l,r'’t-l.3^'’J.3>l"'’t.3-l’ 


[A. 4] 
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The minimization equation (A. 3] is solved by simul- 
taneous over-relaxation. The matrices S . and G. . mav 

be computed initially except for the A. . term and will 

^ » 3 

not change throughout the analysis. Matrix H. . must be 

^ » 3 

recomputed for every iteration of the relaxation. 

The relaxation proceeds as follows; At Point (i,j) 
the terms of the minimization equation are evaluated using 
the assembled P field for P*. In general, the equation is 
not satisfied and a residual is defined as 


S 


1/3 


P* 

1/3 



+ 



R 


[A. 5] 


The superscript t is an iteration counter. The value of 

P^ . is to be altered so that on the next iteration, the 
•^ / J 

residual v/ill be zero, provided H. . does not change. Of 

1/3 

course, H. . will change, but if the equation is fairly 
well behaved, repetition of the procedure should lead to 
convergence on the correct solution. 


S. . P* - (G. . + H. .) = 0 

1/3 1/3 ' 1/3 1/3 


Subtracting [A. 6] from [A. 5], 


s. . (p* T - p* = R 

1/3 1/3 1/3 


[A. 6] 


IA.7] 
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and 


_* x + .T 
1/3 


P* T 
1/3 



Convergence can be hastened by increasing the correction 
term in (A. 7] by a factor ALFA. The factor by which it is 
increased is called the over-relaxation coefficient. 

Equation [A. 7] becomes: 

P* = P» ^ - ALFA 5? [A. 8] 

1,3 1,3 Sj^^. 

One iteration consists of making the correction [A. 8] 
at every grid point. Testing has shown the convergence can 
be speeded up and unwanted solution noise decreased if the 
grid points are processed in a circular manner. Therefore, 
the field is scanned in a counter-clockwise circular sweep 
star^-ing at the center and working toward the boundaries. 
Iterations are repeated until the maximum residual is less 
than a specified convergence criterion. The resulting P* 
field is the solution of equation [A . 3 ] . 
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IV. 


RF;-FVArAIATINr. THR DATA WEIGHTS 


At- the end oC each cycle, the weiqht of each report is 
reevaluated . An observation will have its weiqht reduced if 
the report differs from the analysis value on the current 
scan by more than a subjectively determined amount. REVAL 
is the reevaluation parameter and CRIT, the critical value 
at which a report's weiqht is reduced. REVAL is calculated 
as : 


RE At, - ( onwT*opA'rT)' 

and 

OrA'^'T - REEAr**(FP-'' 

where FP - scan numher 

RRFAG - constant, reevaluation factor 
OnWT oriqinal data weiqht 

niF difference between renort and anal v zed 

field value for that location. 


If REVAL is less 
oi'iqiiial weiqht, 
Stwiu. Tf REVAf, 


nwT 


where 


tluin CRIT, the observat ion retains its 
ovi’u tluHiuh it lUiiy have been reduced another 
is qi'eater tluin CRIT, then: 

con; IT * onwT 
"'l, + REVAL 

CONRT - constant 

'CPWT = o r i q i na 1 da t a we I q h t 

- new d>ita wi'iqht 
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5 - • 'fx-t 

icy? i. 


Notice that on any cycle, every data point may have its 
original weight restored, even if it had been reduced pre- 
viously. In this way, a report that causes a large change 
in the analysis may have full effect if it is supported by 
data nearby. 
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APPENDIX B 


WIND ANALYSIS USING THE PATTERN CONSERVING TECHNIQUE 
I . INTRODUCTION 

The essential feature of the pattern-conserving tech- 
nique is that, while fitting new data, it tends to retain 
certain differential properties of the first-guess field. 
Some of the properties of the wind field we would like 
to conserve; e.g., vorticity and divergence, involve both 
scalar wind components. 

The differential properties that we choose to conserve 
are the gradients of each wind component in eight directions 
from each grid point, the vorticity and the divergence. 

The same method is used here as in the scalar analysis, 
the main difference being that two minimization equations 
rather than one must be solved simultaneously. 

The equations are simplified by using the staggered 
grid illustrated by Figure B-1 and defining the divergence, 
vorticity and gradients as in Table B-1 and Figure B-2. 

This arrangement causes certain matrices to be tridiagonal. 
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II. ASSEMBLY 


The method of assembling data to grid points is similar 
to that in the scalar analyses. For the wind analysis, 
the u and v components are treated separately except for 
the gross reject criterion which is a function of the 
guess-field wind speed. The assembly radius for a particular 
report is based upon the cycle number and the information 
density at its location (as in the scalar analysis) . The 
function which determines the weight of an observation is 
inversely proportional to both the local vorticity and to 
the relative distance between the observation and the grid 
point. Note that since the grid is staggered (Figure B-l) , 
the individual wind components of an observation may be 
assembled to grid points with different weights. 
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III. MINIMIZING THE DEVIATIONS 

The objective of the Pattern-Conserving Technique 
is to preserve the differential properties of the guess 
field. In the scalar analyses, it was found that to use 
the finite difference expressions centered about one grid 
point did not provide enough horizontal coupling. Accord- 
ingly, the integral was expanded to include the expressions 
for the surrounding grid points as well. In the wind analy- 
ses, this is not necessary because the particular dif- 
ferential properties and the staggered grid extend the 
influence of the observations more than in 
analyses and facilitate coupling. 


the scalar 


TABLE B-l 


PCT VECTOR CONSTRAINTS 


Constraint 

u. = Variable being analyzed (assembled value) 
1 


V feight 

^l,m 


^ = Variable being analyzed (assembled value) 


‘l,m 


l,m 


divergence = + 3v/ 

o X 


3y 


'■^ 1+1 ,m “ ^l,m ^l,m+ i '^Ijm 


*1 ,m 


(Computed from non-assembled value of first guess.) 
vorticity = - 3u/g^ 


r . 

1 ,m 


1 ,m 


1 ,m 
®1, m 
^l,m 
^l,m 
^ 1 , m 
*^l,m 
^l,m 
^l,m 


V. - V, , - u, + U-, 

l,m 1- i,m l,m l,m-i 

(Computed from non-assembled value of 

-x-i,v+i u gradient = u, , - u, 

^ 1- 1 ,m+ t l,m 

(Computed from non-assembled value of 

X- 1 y+ ) V gradient = v, , , , - v-. 

’ 1- 1 i,m 

(Computed from non-assembled value of 

y axis u gradient = u, , ~ u, 

^ l,m+ 1 l,m 

(Computed from non-assembled value of 

y axis v gradient = v, , , ~ v. 

^ l,m+ ‘ i,m 

(Computed from non-assembled value of 

x+i,y+i u gradient = u,. , , , - u, 

■' ^ 1+1 ,m+ 1 l,m 

(Computed from non-assembled value of 

x+i,y+i V gradient = v, , . - v, 

^ ^ 1+1, m+1 l,m 

(Computed from non-assembled value of 

x axis u gradient = u ~ ^ 

1+ I ,m l,m 

(Computed from non-assembled value of 

X axis V gradient = v, , - u, 

^ 1 + 1 ,iTi l,m 

(Comouted from non-assembled value of 


first guess.) 
first guess . ) 
first guess.) 
first guess . ) 

first guess.) 
first guess.) 
first guess . ) 


.,m 


' l,m 


l,m 


1, m 


l,m 


■L,m 


H 


l,m 


first guess . ) 

H 

first guess . ) 


l,m 
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We shall minimize the following integral; 


t = I! t 


(a) 


A, „(u* „ - u, ) 
i,in 


I H.IJ 


(b) 

(c) 


+ A, (v* - V, ) 




(d) 


- n-i.r. ~ '^I,rn 


(e) 

+ 


,m ^'^ 1 - i ,m+ 1 '^1 







(f) 

+ 


( V * • V* 

,m 1 - i ,m + 1 1 ,m 

■ 

(g) 

+ 


(u* - u* 

,m 1 ,ir 1 1 ,n\ 


(h) 


A 


A 

+ 


(V* ^ • V? 

,m 1 ,m + 1 1 ,m 


(i) 

- 1 - 


,n ^^ 1 + 1 ,m+ 1 ,m 

- 

(J) 

+ 


(vf - vf 

,n' l+i„m+i l,ni 

• * 31 ,, 


(k) 

( l ) 


+ H, (u* "U? ”h,„) 

1 ,m 1+ I ,m 1 ,m 1 ,r» 


+ H, (v* - V* - h, ) 

1 ,m 1+ 1 ,m 1 ,m 1 ,m 


J dxdy 


The superscript (*) indicates the values we seek. The 
differential properties of the first-guess field are defined 
in Table h-l, and a depiction of the u component minimization 
stencil as it relates to the u terms of equation I B .1] is 
given in Figure B -2. 

To minimize the integral we take the first variation with 
respect to uj ^ and with respect to v* yielding the 
following two equations: 


61 

•//I 

A, (uf " u, ) 

l,m I.m l,m 


( B -21 



,m 

'"I*,, 

- u? * V? ” V? 

1 . .T* 1 , m+ j 1 , m 

- "l,.’ 



1 , z\ 




- 


n+i " ^l,n ■ 


- 



1 ” ” ^1, m 



_6j 

6v? 

1 ,m 

= // 

^ ^ 1 , m , m '^1 ,m^ 


1 B .3] 

^l,m 


-u* +v? ^ -V* 

l,m l,m+j 

- ‘'l.m’ 


^ , m 


-V*. "VI? +u? 

i- 1 , m 1 , m i , m - 1 

- ^1,:.’ 



- 

(v*i_ 

- }/* - 
i > m + 1 1 , m 

e, ) 
1 ,m 

" ^l,m 

(v? - V* - f , ) 

1 . m-t 1 1 , p. 1 , m 

- 


- V ? - 

, m+ 1 1, m 


- »1... 

i+i,m l,p l,rti = 
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In equation [ B .2] group terms involving 1) uj 2} u* at 

surrounding points; 3) v* and 4) everything else. 


’l.m 


t B.4) 


fn(A, „ ♦ D, + Q, + E, + P, ♦ G, ♦ H, ) U* 
ji l.m l.m l.m l.m l.m l.m l.m' l.m 


U.m’ 


* ’■ °l.n ■ '"Ol.a’ “l.n-, 


. m 


Y, (D, - Q, ) V* + (-D, ) V? ♦ Q, V? 

l.m l.m '^l.m i.m ' l.m’ l.m-*-i '^l.m ^l-i.; 

^l.m '^l.m ’^l.m '^l.m °l.m '^l.m ^1 . n ® 1 , m * ^1 .n ^l.m 

* ^l.m -^l.m ^ »l.m ^l.m’ = ° 


'l.m 


Group [ B .3] similarly: 


l.m 


I .51 


//[{A, „ + D, „ + Q, +E, +F, G, -*• H, ) V* 

'' 1." l.m l,m l.m l.m l.m l.m 1, 


m 


‘1 . m 


r+ (-D, - F. V* +(-Q, „ + (-E, ) V* 

l.m 1**^ l.m+; l.m l-i , n l.m 1-i.m+i 


^l.m ^^l.m ^l.m^ '^i.m ^l.m '^l+>.m "*’*^l.m ’^l.m-i 


‘l.m 


T" ^l.m ''l.m'^^l.n ‘^l.rn'^l.m ‘^l.m'^ ^1 .m® 1 . m ■*■ ^1 .m ^l.m 
^ ^l.m ^l.m-^»l.m ^l.m’ ° 


Note that all terms in S and Z except A, in S, and 

l.m 1 , m 

- A, u- in Z, involve first-guess information which is 
x^rnl^ru X^ri 

consta"*t during tr.e . r.alysis. Similar conditions hold for 

A A 

S and Z . 
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Equations ( n.4] and t ^.5] can be written in matrix form: 


“1 ,m 

u* 

+ X, 

•“i ,m 

^ il,m 

"*■ 2. = 0 
"1 ,m 

1 B .6] 

A 

“1 , m 

V* 

A 

+ X, 

A 

+ Y, 

+ Z, =0 

r B .7] 


-“1 ,m 

~1 ,m 

-1 ,m 


These equations must be solved simultaneously. The method 
of solution used is Liebmann successive over- relaxat ion. 
Using a first-guess for u* and v*, equation [R.6] isr in 
general, not satisfied, A residual is defined by: 


S, + X 

•“i ,m — “• 


, + Y, + Z, = R 

1 ,m “■! ,m *^‘1 ,m 


[ B .8] 


The superscript t is the iteration counter. We wish to find a 
next guess at u* such that the residual is zero, if the values 
at surrounding points do not change. 


S, + X. + Y. + Z, = 0 

“1 ,m— “I,::: “l,m *=l,m 


[ B .9] 


Subtracting [ 0 .9] from ( 0 .8] , 


£l,„ (i*' - = K 


and u*''-'' .. U*'' - 


> , 

=1 ,m 


I B.IO] 
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Convergence is more rapid if the correction in I W .10] is 
augmented by the inclusion of ALFA factor. 


u*'^'*"* = - ALFA . 

-l,m 


I B .11] 


At a particular grid point, u* is corrected by 
equation [B.ll] and v* is then corrected in an analogous 
way. In computing R from equation [ B .8] or from the 
analogous equation in v* , the latest estimate of both u* 
and V* at surrounding points is used. Some of them have 
been changed on the current iteration and some have not. 

As in the scalar analysis, the field is scanned in a 
counter-clockwise circular sweep starting at the center 
and working toward the boundaries. 

During each iteration through the grid, the maximum 
residual is ch^cVed. When it becomes less than a pre- 
scribed convergence criterion, the equations are considered 


solved . 




IV. RE-EVALUATING THE DATA WEIGHTS 

The validity of wind reports is judged according to the 
vector difference between the reported wind and the analyzed 
wind. The analyzed wind is obtained by interpolation from 
the analysis fields. Tne reevaluation parameter (REVAL) 
is the tame as its counterpart in the scalar analysis 
with the exception that DIF is the magnitude of the vector 
difference squared, i.e. 

DIP = IWj-W^P 

where \V^ = interpolated analyzed wind 
= wind report 

If REVAL is greater than a specified critical value, the 
weight is reevaluated in the same way as for the scalar 
analyses (see Section IV of Appendix A) . Again, if REVAL 
is less than the critical value, a report is assigned its 
original weight, even if it has been reduced on a previous 


scan. 


APPENDIX C 


^^ASS STRUCTURE LINEAR TRANSFORMATIONS 


I. INTRODUCTION 

If an analysis of upper-air pressure heights and tem- 
peratures is to be used in initializing a forecast model, 
it is desirable for the heights and temperatures at grid 
points to be consistent with the hydrostatic equation. It 
will be shown in Section II that the heights, temperatures 
and layer stabilities can be interrelated through various 
linear transforms. It turns out that to close the equation 
set it is also necessary to specify a single height parameter 
and a single temperature parameter. 

The vertical organization of height and temperature level 
and stability layers to be used in the mass structure analysis 
is shown in Figure C-l. The stability parameter used here is 
defined as: 

2 

p 31n9 r « 1 1 

c I C. 1 J 

P 3P 

Other definitions are possible, as discussed by Holl ^ ^ 
(1963) . This definition makes 0 linear in p i<= R/Cp) , 
which is consistent with pseudo-adiabatic diagrams. 

A limitation of this technique is that o must be assumed 
to be c instant in each of the layers labeled 1 - 10 in Figure 
C-l. If a serious departx’re from this condition occurs in 
a layer, the temperature above the layer will depart from 
the reported temperature, but will agree hydrostatically 
with the analyzed heights. 



II. DERIVATION OP THE TRANSFORMATIONS 


The hydrostatic 

dz _ 1 

dp ” gp 


equation states 


RT 

pg 


where p 



or 

T = - E£ ^ 

R dp’ 

Potential temperature is defined as 


P 

G H T (p-) where < 
o 

Therefore , 



G 

Def ini 
In 0 = 


qdz l-< 

R dp P <Po' 


ng n - l-<, 
in I t In (-||) 


Ino +<lnP 

o 


I C.2 ] 


I C.3 ] 
( C.4 ] 


[ C.5 ] 
I C.fi 1 


and 

= d Hn(-^) ]+ a I C.7 ) 

dp dp dp p 

and 

2 

dz dlnO _ d z H dz r r> o i 

•• "p j ''' ”* *“ ■ • I ^ ^ j 

dp dp p dp 

Substituting for p from the hydrostatic equation 
into [ C. 1 ] gives 


a 



dlnO 

dp 


I C.Q 1 


• C-3 


and 


dln0 


dz 2 • 


I C.IO ] 


Substituting { C.IO ] ..nto { r'.fi ] 


cl z ^ ^ 

P dp 


Multiplying by 


n d z 


+ n? 


n-i dz 


2. pH-n 
g ^ 


Let us define 


X H - p 


0 

dp 


Then 


^ = 
dp 




From ( C.12 ] and ( C.14 ] 


- ^ - o „n-2, 

dp - 5 P 


Integrating within the layer characterized by 
constant stability 


- X + c, . 2 p''-^ + 

1 g n~l 


t C.16 ] 


Defining M = C2 - 


X = p'’"^ + “• 


Integrating [ 0.1“^ 1 
Z ® C-, - / dp . 


Substituting ( C.17 ] into ' C.18 ) 


-<=3-^ f >57^1 P‘" - "P'"! <5P 


^g(n-i) 


or 


z =c,+ Inp - M ^ 


■3 g(n-l) 

and finally 


Z = N* + M* p*^ - — In p 

g< ^ 


where 


N* = 


and 


M* = - 


M 


1-n • 


( c.iM 


[ C.18 ] 


[ C.19 ] 
( C,20 ] 


[ C.21 ] 


Equation [ C. 21 ] is the basic equation of the method 
relating pressure heiqht to stability . We also need a 
relationship between the temperature arid the stability. 
Taking the first derivative of equation [ C.21 J with 
respect to pressure. 


‘^2 . K-l o 

- M* K p - . 

dp ^ g*^P 


I C.22 ] 
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Substituting [ C.22 1 into I C. •? ] 


T = - |SL (M* ►: 




[ C.23 ] 


or 



M* I- p^ 


I C.24 ) 


Equation I C.24 ] is the basic equation of the method relatinq 
temperature to stability . 

Equations [ C.21 ] and ( C.24 ] are the two model equa- 
tions we need. They apply to each of the ten layers in Figure 
C-1. The N*, M* and a in the ten layers make a total of 
30 unknowns. Applying equation ( c.2l ) to each mandatory 
level gives us twelve equations: 


”1 + «1 '■a - “l®a = 

”1 * "l '’b" - 


1 C.26 1 


N + K P, - o B, = Z, 

10 10 k 10 k k 

N, + M, P, - a, S, = Z. 

10 10 1 10 1 1 


where 


3 = — In p and the subscript (*) has been omitted, 
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1 

I 



Requiring continuity of height at the interface 
levels between each layer leads to nine more equations: 


V '"l - ”2> - '"l 




I C.27 J 


»9 - N„ + Py' (Mj - M„1 




- 0 ») = 0 


Requiring continuity of temperature at the interface 
levels gives, from equation I C.‘?4 ], the remaining nine 
equations : 


“m '”l - «2> - '“1 - »2> ■= 0 


1 C.28 ) 


Oy (Mg - Mjj) - (Cg - Qjg) 


= 0 


where 

2 ic 

a = g< p , 

The 30 equations in 30 unknowns may be written as a 
single matrix equation. 

B C = Z*. I C.29 ] 


The vector is composed of the twelve mandatory level heights 

and 18 zeroes. The vector C is the 30-element column vector 
N 

(M) , where the ten elements of N, and M and ^ correspond to 
a 
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the ten layers. Equation [ C.29 ) is written out in Figure 
C-2. This can be represented in a partitioned form as 


^1 ^2 ^3 


N 


z 



A 



Hi H3 


M 

S 

0 



A 



Gi G2 C3 


0 


0 


with the formal solution 
C = Z' . 


v/hich can be represented in partitioned form as 



[ C.30 J 


( C .31 ] 


In the analysis/ we need a transformation to get stabi- 
lities from heights. That transformation is part of matrix b" 
namely 2 . = Z. 

We will also need a transform back to heights. For that 
problem, our set of 30 equations contains 32 unknowns (10 Ns, 
10 Ms and 12 Zs) . Two of the unknowns will have to be given 
in order to close the set. 

The obvious choice for one of them is the 1000 mb height, 
since more data is available at the surface than in the upper 
air. Choosing the second parameter is more difficult. Since 
the temperature will be computed from the heights, the second 
parameter might be chosen as a reference for the temperature 
profile. We chose the thickness of the 1000 - 250 mb layer. 


0 10 
0 0 1 
0 0 0 


0 0 0 


0 0 0 


0 1-1 


0 0 0 
0 0 0 


0 0 0 
0 0 0 


0 0 0 
0 0 0 


0 0 0 


0 0 f, 0 0 P^O 0 0 0 0 0 0 0 0-0^0 00000000 
0 0 000 p !^0 0000000 0 - 6 . 0 00000000 
000000 pIo 000000 0|0-S 00000000 


0000 0.0 0 PjO 0 0 0 0 0 0|0 0-0^0 000000 

0 0 0 0 olo 0 0 p‘‘0 0 0 0 0 0.0 0 0-6^0 0 0 0 0 0 


0 0 0 0 OlO 0 0 0 PjO 0 0 0 OlO 0 0 O-0J^O 0 0 0 0 

1000000000 0*^0 0 0 olo 0 0 0 c-e 0 0 0 0 

9 , ' 9 

0 10) cjo 00000 pJ^O 0000000 0-0^0 0 0 
0010 ojo 0 0 0 0 0 0 pfo ojo 00000 0 - 3.0 0 
0001000000000 pI'oIo 0 0 0 0 0 0 O-B^O 

0 0 0 0 110 00000000 p|;o o o o o o o o o-$. 


0 o|o 0 0 0 0 0 0 p:o olo o o o o o 0-3.0 0 


000 0 1|0 OOP 0 000 0 PjO 0000 O OP 0 -S, 

0 0 0 0 0 P“'’-P'‘O 0 0 0 0 0 0 0-3 BOOOOOOOO 
m in ^ mm 

0 0 0 0 0 0 P*^-p “'0 0 0 0 0 0 010-6 80000000 
n n Inn 

0 0 0 0 0 0 0 P‘‘-P 0 0 0 0 0 olo 0 -B 8 0 0 0 0 0 0 
0 ^ 0 ^ I ® O 

0000000 OP -P 0 0000000-0000000 

p p I p 

-1 0 0 0 0,0 0 0 0 P‘‘-P *"0 0 0 o |0 0 0 0-3 3 0 0 0 0 
I ^ ^ I s ^ 

1-1 0 0 ojo 0 0 0 OP^-P^O 0 ojo 0 0 0 0 - 6 j.Sj .0 0 0 

0 1-1 0 o'o 0 0 0 0 op'^-p'^o 0,0 0 0 0 0 0-3 6 0 0 

s .s I s s 

0 0 1-1 OjO 0 0 0 0 0 op^-p^olo 0 0 0 0 0 O- 3 ^ 0 J.O 

0 0 0 l-l|o 0 0 0 0 0 0 cp|^-p^j 0 0 0 0 0 0 0 0 - 8 jjS, 

000 0 Oa -u 0 0 000 ' 0 0-1 1000 000 3 0 


OOa-uOOOOO'O 0-1 100000000 
m m 

0 OiOo -a 0 0 0 0 0 0 010-1 10000000 

Inn I 


0000 OjOo^-a^O 00000 010-1 10000000 

0000 ojo Oa^-a^O 0 0 0 0 ojo 0-1 1000000 
0000 o'o 0 Oa -o 0 0 0 0 0*0 0 0-1 100000 

p p I 

0000 010 0 0 Oa -a 0 0 0 010 0 0 0-1 10000 


0 0 0 


0 0 0 
0 0 0 


0 ojo 000 Oa^-n^O 0 0|0 000 0-1 1000 

0000000 Oo-ct 00, 0 0000 0-1 1 0 0 

■ s s I 


0 ojo 0 0 0 0 0 0aj,-aj,0j0 0 0 0 0 0 0-1 1 0 
00 lo 000000 Oa^-a^O 0 0 0 0 0 0 0-1 1 


FIGURE C-2: 





Oefino the 12-element column vector E = I H 


\ 1 / 

where Z is the iOOO mb height and ii is the 1000 - 250 mb 
thickness. We need the trans 'ormation E = D Z. 

The last ten rows of matrix D are the first twelve 

8S 

columns of the last ten rows of of ( C.31 3). The 

first two rows of D are 


100000000000 

-lOCOOOOOlOOO. 

Matrix D ^ may be obtained by Gauss elimination, and the 
heights can be recovered using ^ £. 

Let us repeat the sequence of operations. At grid 
points, matrix g.^ is used to convert twelve mandatory level 
heights to ten layer stabilities. The stabilities are limited to 
be greater than zero and less than a maximum value if necessary. 
Then matrix D ^ is used to compute the mandatory level heights 
at the grid points. 

The temperatures at the grid points can now be computed 

by simply substituting o and M* , which are submatrices of 

** X ^ 

B , into I C.24 ]. In matrix form, T = Q Z where T com- 

A 

prises the twelve mandatory level temperatures, Z the twelve 
mandatory level heights and Q the matrix of [c.24] 
namely 

Dl/R< - 2- E^. 

P 


[ C.32 1 



The ter.peratures and heights at the layer interface 
levels car be obtained by simply changing the pressures 
in the matrix coefficients. By proper substitution, the 
following matrix equation can be formulated: 


where 




Vectors and Zj are the nine temperatures and heights at 
the layer interfaces. 
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APPENDIX D 


RADIOSONDE CHECKER 

Radiosonde soundings and satellite temperature/height 
profiles are known to occasionally contain errors due to 
a garbled transmission or errors in the actual workup of 
the sounding. The hydrostatic equation can be used as 
a check on whether reported values at successive levels 
in the sounding are reasonable. If levels are missing 
in the sounding, the hydrostatic equation can be used to 
approximate a value for them. In the current model the 
heights and temperatures at the 950 and 900 mb levels are 
analyxed. However, these are not mandatory levels for 
radiosonde reports, so all ’’observations'' must be generated 
via this interpolation scheme. 

The first part of the radiosonde checking program 
builds the working array. Significant , and then mandatory , 
level observations are read and duplicates are removed. 

The mandatory and significant level observations are then 
merged into an array in which pressure monotonically de- 
creases (height increases) . Tropopause and maximum wind 
data, if reported, are inserted at the appropriate level. 
Three major steps complete the process; the array is checked 
for hydrostatic consistency; heights are assigned to signifi- 
cant level reports; and missing levels are approximated 
(if possible) . 




In the first part of the hydrostatic consistency evalua- 
tion, the temperature lapse rate is checked using 3 degrees 
C per 100 meters as a gross measure. If the lapse rate 
at any ]evel exceeds this, the level is flagged as missing 
and the check continues until the top of the sounding is 
reached. Next, an attempt is made to find the station 
level report to extract the pressure information. Generally, 
this is the first level reported in a significant level 
report. Knowing the station elevation, one can calculate 
the standard station pressure. This is compared to the 
reported station pressure and, if the two values differ 
significantly, it is assumed that the station level report 
is not available. If this is the case, the first reported 
mandatory level is used for the base height for the hydro- 
static workup. 

Given a surface pressure, height and temperature, 
one car calculate, using the hydrostatic equation, the 
height at the next pressure level if the pressure and temper- 
ature at that level are known. These heights are included 
in the significant level reports and changed in mandatory 
levels if the report does not appear to be consistent. 

Below 250 mb, the average height char.ge required to make 
the sounding hydrostatically consistent is about 5 meters 
vln a sample size of approximately 450 soundings). Above 
250 mb, the change varies from about 8 meters to 10 meters. 
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Next, missing information in the soundings is replaced 
by interpolated values. This includes heights, tempera- 
tures, dewpoint depression, wind direction and wind speed. 
The missing levels are interpolated in ln(p) wherever pos- 
sible. V7ind reports from both mandatory and significant 
levels are merged with tropopause and maximum wind reports. 

Finally, the sounding data for the levels analyzed 
by the mass structure programs are extracted and written 
to the disk in a format compatible with the analysis input 
requirements. 
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STRATOSPHFRir HKir,HT-TE«PFRA’’’URF EXTRAPOLATION 

The method for obtaining stratospheric height and 
temperature fields is based on the work of Lea (ThfiT). 

The extrapolated values are given bv equations of the form 

==V*l^level-l*^2'^level-l 

In this manner, the 'SO mb fields are extrapolated from 
the 100 mb height and temperature, 10 mb from the *S0 mb, 
and the 10 mb from the 10 mb. The coefficients, which 
are a function of latitude and month, were obtained from 
empirical studies of selected rawinsonde stations. After 
the fields have been produced they are filtered. 
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